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AFIT/GAE/ENY/97D-05 


Abstract 


Recent  research  efforts  have  applied  the  receding  horizon  Model  Predictive  Control  (MPC) 
strategy  to  linearized  high  performance  aerospace  systems.  The  research  contained  in  this  the¬ 
sis  used  these  recent  results  in  order  to  apply  the  MPC  strategy  to  a  nonlinear  high  performance 
aerospace  system,  specifically  an  F-16  fighter  aircraft  model.  The  model  was  commanded  to  fol¬ 
low  dynamic  trajectories  of  roll  angle  and  altitude.  Further,  adaptive  constraint  techniques  were 
used  to  improve  system  tracking. 

To  accomplish  these  tasks,  code  and  block  diagrams  were  generated  using  the  commercial 
software  packages  of  Matlab  and  Simulink.  Numerous  simulations  were  conducted  with  the  goal 
of  achieving  realistic  aircraft  performance.  In  many  cases,  to  improve  system  tracking  and  reduce 
control  input  oscillations,  rigid  mathematical  constraints  previously  used  in  the  MPC  strategy  were 
relaxed. 
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CONSTRAINED  MODEL  PREDICTIVE  CONTROL 
OF  A  NONLINEAR  AEROSPACE  SYSTEM 


Chapter  1  -  Introduction 

1.1  Model  Predictive  Control 

Model  Predictive  Control  (MPC)  is  a  control  strategy  which  optimizes  a  specified  performance 
index  over  a  set  of  future  inputs  to  minimize  future  output  deviations  from  a  specified  trajectory. 
MPC  typically  operates  on  three  receding  horizons:  an  optimization  horizon,  a  control  horizon  and 
a  prediction  horizon.  An  on-line  optimization  takes  place  over  the  optimization  horizon  while  future 
control  inputs  are  calculated  over  the  control  horizon  and  future  output  trajectories  are  calculated 
over  the  prediction  horizon.  The  controller  then  implements  the  first  control  input,  discards  the  rest, 
and  recalculates  the  next  series  of  control  inputs  at  the  next  time  step.  Because  of  the  computational 
intensity  due  to  the  on-line  optimization,  MPC  has  traditionally  been  applied  to  low-bandwidth 
processes. 

Recent  work  has  focused  on  applying  the  MPC  strategy  to  high  performance  linearized  aerospace 
systems  [4] ,  [11] ,  [6] .  Frequently,  MPC  formulations  which  use  a  stabilizing  inner  feedback  loop 
have  been  applied  to  such  systems.  The  use  of  a  stabilizing  inner  feedback  loop  guarantees  system 
tracking  within  a  finite  number  of  time  steps  for  discrete  controller  poles  placed  at  the  origin.  Future 
commanded  inputs  are  generated  by  optimizing  over  a  reference  signal,  v,  which  is  the  input  to  the 
inner  feedback  loop.  This  optimization  allows  for  an  independent  optimization  horizon,  r,  because 
the  optimal  plant  inputs  are  no  longer  directly  calculated.  Also,  the  use  of  the  inner  stabilizing  feed¬ 
back  loop  ensures  stability  through  the  existence  of  a  monotonically  decreasing  cost  function.  The 
distinct  advantage  of  the  MPC  strategy  over  traditional  controller  designs  is  its  ability  to  account  for 
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changing  system  constraints,  i.e.  actuator  limits  and  rates,  as  well  as  output  state  constraints  such 
as  altitude  and  bank  angle. 

1.2  Importance  of  Research 

Because  the  MFC  strategy  is  capable  of  dealing  with  changing  constraints,  it  is  potentially 
suited  to  aerospace  systems  in  which  complete  or  partial  actuator  failures  are  possible.  With  aircraft 
which  are  statically  unstable  or  at  best  neutrally  stable,  human  response  may  not  be  sufficient  to 
maintain  positive  control  of  an  aircraft  following  an  actuator  failure.  Also,  as  flight  research  con¬ 
tinues  into  the  use  of  new  combinations  of  actuators  and  reduced  torsional  stiffness  wings  [10]  , 
designing  controllers  which  provide  the  greatest  control  authority  using  the  least  control  power  may 
be  difficult  using  classical  control  techniques.  The  MFC  strategy  is  also  suited  to  flight  test  work 
because  of  its  ability  to  handle  changing  constraints.  In  an  actual  flight  test,  the  MFC  controller 
could  potentially  act  as  safety  watch  dog,  preventing  the  pilot  or  test  from  exceeding  certain  para¬ 
meters.  By  predicting  these  future  parameters,  the  actuator  constraints  can  be  modified  real  time  to 
prevent  a  dangerous  situation  before  it  occurs.  MFC  can  also  be  coupled  with  current  Filot  Induced 
Oscillations  (FIO)  research.  By  predicting  future  outputs,  pilot  inputs  can  be  modified  to  assist  cur¬ 
rent  filter  designs  in  the  prevention  of  FIO. 

1.3  Research  Objectives 

Before  MFC  can  be  realized  on  an  actual  aerospace  system,  it  must  be  shown  to  work  well  on 
simulations  involving  a  nonlinear  plant.  Building  upon  past  work,  this  thesis  examined  four  areas  of 
use  of  an  MFC  controller.  The  first  was  to  expand  the  recent  work  [4]  of  step  responses  and  demon¬ 
strate  MFC’s  ability  to  handle  dynamic  trajectories.  The  second  was  the  application  of  the  MFC 
strategy  to  a  nonlinear  F-16  aircraft  model.  The  third  area  was  the  relaxation  of  rigid  mathemati¬ 
cal  constraints  associated  with  pole  placement.  And  last  was  the  application  of  adaptive  constraints 
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to  improve  system  tracking  and  reduce  overall  required  control  power.  These  four  areas  were  re¬ 
searched  with  an  overriding  emphasis  that  simulations  should  model  realistic  aircraft  performance. 

1.4  Thesis  Overview 

Chapter  2  explains  the  mathematical  derivations  of  the  state  space  MFC  formulation.  Devia¬ 
tions  from  earlier  works  are  noted,  and  potential  tuning  parameters  are  identified.  The  chapter  also 
addresses  the  rigid  mathematical  constraints  previously  used  in  MFC  studies  and  offers  a  heuristic 
approach  to  achieve  asymptotic  stability.  Chapter  3  details  the  F-16  linear  and  nonlinear  aircraft 
models  which  were  used  in  this  research.  The  two  models  as  well  as  the  entire  MFC  block  dia¬ 
gram  are  described  in  detail.  Additionally  flight  conditions,  open  loop  pole  locations  and  potential 
problems  associated  with  differences  between  the  two  models  are  identified.  Chapter  4  presents  the 
results  of  various  Matlab  simulations  which  were  designed  to  demonstrate  MFC’s  ability  to  meet 
the  research  objectives.  Chapter  5  offers  conclusions  and  recommendations,  based  upon  the  results 
of  the  research,  providing  starting  points  for  future  research. 
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Chapter  2  -  Review  of  Literature  and  Mathematical  Development 


2.1  General  Predictive  Control  and  Stable  Generalized  Predictive  Control 

This  section  covers  the  basics  of  two  of  the  more  common  types  of  Model  Predictive  Control. 
The  first  is  General  Predictive  Control  (GPC)  and  the  second  is  Stable  Generalized  Predictive  Con¬ 
trol  (SGPC).  The  advantages  and  disadvantages  of  each  are  discussed  and  a  suitable  method  was 
chosen. 

2.1.1  Generalized  Predictive  Control 

GPC’s  basic  concept  is  to  minimize  a  cost  function  across  separate  input  and  output  horizons, 
where  a  horizon  is  an  integer  number  of  steps  into  the  future.  Over  the  output  horizon,  the  tracking 
error, 

e  =  y{k)-s  (k) 

y{k)  =  System  Outputs  at  Hme  Step  k  (1) 

s{k)  =  Commanded  Trajectory  at  lime  Step  k, 
is  minimized.  Over  the  input  horizon  the  optimization  minimizes  the  required  control  power  to 
reduce  the  tracking  error,  e.  To  do  this,  GPC  uses  an  explicit  plant  model  of  the  form 

x{k  +  1)  =  Ax{k)  -f  Bu{k) 

y{k)  =  Cx{k) 

to  predict  the  plant  outputs,  y,  across  the  output  horizon,  N.  The  control  inputs  are  weighted  and 
optimized  over  a  separate  input  horizon,  iV„.  The  Single  Input  Single  Output  (SISO)  expectation 
cost  function  is 

{N  Nu 

Y^[y{k  +  l)-s{k  +  l)]‘^  +  xY^  Au{k  +  1-1)2 
1^1  1=1 

Au{k)  =  u{k)  —  u{k  —  1) 


(2) 


where  A  is  the  weighting  of  the  control  power.  An  expectation  operator  is  used  because  the  output, 
y{k  + 1),  can  be  contaminated  by  various  noise  sources. 

There  are  several  advantages  of  GPC  over  other  controllers.  First  there  is  no  requirement 
to  know  the  system  closed  loop  poles  and  there  are  no  ill  effects  caused  by  closely  spaced  zeros 
and  poles  [4] .  GPC  also  allows  the  use  of  control  input  constraints  as  well  as  system  constraints. 
These  constraints  provide  the  control  designer  more  tuning  parameters  to  improve  system  tracking. 
However,  GPC’s  major  fault  is  it  offers  no  guarantee  of  stability. 

2.1.2  Stable  Generalized  Predictive  Control 

Stable  Generalized  Predictive  Control  (SGPC)  builds  upon  GPC’s  advantages  and  adds  a  stabil¬ 
ity  guarantee.  This  stability  guarantee  is  realized  by  the  introduction  of  a  stabilizing  inner  feedback 
loop  based  upon  the  Youla  parameterization  of  all  stabilizing  controllers  [7]  .With  the  introduction 
of  an  inner  feedback  loop,  the  MPC  controller  now  operates  on  a  reference  input  v{k).  This  allows 
the  introduction  of  a  new  optimization  horizon  r  over  which  the  cost  function  is  minimized.  If  the 
discrete  poles  of  the  stabilizing  inner  feedback  loop  are  chosen  to  be  at  the  origin.  Finite  Impulse 
Response  (FIR)  behavior  is  achieved.  FIR  behavior  implies  that  a  system  will  achieve  a  steady  state 
value  over  a  finite  horizon.  When  used  with  a  cost  function  that  minimizes  tracking  error,  both  sta¬ 
bility  and  a  monotonically  decreasing  cost  function  are  guaranteed  provided  specific  horizon  con¬ 
ditions  are  met  [4]  .  In  this  thesis  the  requirement  for  FIR  behavior  was  relaxed  and  it  was  shown 
that  a  finite  error  could  be  reached  over  a  finite  horizon.  Because  of  the  ability  to  reach  a  finite  error 
and  the  stability  guarantee,  the  SGPC  method  was  chosen. 

2.2  State  Space  Formulation  of  MPC 

The  state  space  formulation  of  SGPC  was  done  in  a  four  step  process.  First  the  SISO  cost 
function  of  Equation  2  was  rewritten  in  multiple  input  multiple  output  (MIMO)  form.  Once  the  cost 
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function  was  identified,  step  two  involved  developing  a  stabilizing  inner  feedback  loop.  In  step 
three  the  state  prediction  matrices  were  developed,  which  were  used  by  the  optimization  routine  to 
predict  future  outputs,  states  and  inputs.  Finally,  a  Matlab  quadratic  optimization  routine  was 
applied  to  the  equations  developed. 

2.2.1  Cost  Function 

An  expansion  of  the  SISO  cost  function  of  Equation  2  in  MIMO  form  gives 

m  =  E  + 0  -  »{* + i)  III. + E II  A<.(fc + i  - 1)  111.  (3) 


Ry 

> 

0 

Ru 

> 

0 

(4) 

Ry 

= 

Ry 

Ru 

= 

Ru- 

The  selection  of  Ra  to  be  positive  semidefinite  implies  that  infinite  control  power  could  be  com¬ 
manded  of  the  system.  However,  this  was  avoided  as  the  control  inputs  were  constrained.  Due  to 
the  stabilizing  inner  feedback  loop  a  reference  input  v{k)  is  introduced.  This  allows  the  cost  func¬ 
tion,  J{k),  to  be  minimized  over  the  optimization  horizon  r  at  each  time  step  k.  The  difference 

i 

of  system  outputs,  y,  and  setpoint  trajectories,  s,  are  weighted  by  iiy  over  the  prediction  horizon, 

p.  Incremental  control  inputs,  Au,  are  varied  along  the  control  horizon,  q,  while  weighted  by  the 

1 

matrix  RZ  ■ 

2.2.2  Stabilization  of  a  ETI  Plant  using  Coprime  Factorization 

The  controller  used  in  this  thesis  was  based  upon  an  inner  loop  feedback  which  provided  sta¬ 
bility  of  the  plant,  G.  Given  that  the  plant  is  a  proper  and  real-rational  transfer  matrix  that  is  both 
stabilizable  and  detectable,  it  can  be  factored  into  a  right  coprime  factorization  (RCF)  and  a  dual 
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left  coprime  factorization  (LCF)  [14,  127] . 


G{z)  =  N{z)M-\z)  RCF 

G{z)  =  M-'^{z)N{z)  LCF 

These  factorizations  satisfy  the  Bezout  Identity  [14,  126]  in  Equation  6. 

fi  -ru  ri=/ 


■  X 

-Y  ■ 

N 

Y  ' 

.  M 

N 

-M 

X 

From  the  factorization  and  the  Bezout  Identity,  the  set  of  all  proper  controllers  which  achieve  internal 
stability  are  parameterized  by  the  following  [14,  323]  : 

K  =  {-X  +  MZr){Y  +  NZr)~'^ 


K  =  (y  +  ZiN'^  ^(-X  +  ZiM)  (7) 

X,Y,X,Y,M,N,M,N,Zr,Zi  G  RH^. 

Before  using  Equations  6  and  7  it  is  necessary  to  first  define  the  discrete  plant  G  and  introduce  a 
feedback  variable  v. 

°  "  [-0\ 

x{k  +  1)  =  Ax{k)  +  Bu{k)  (9) 

y{k)  =  Cx{k)+Du{k)  (10) 

Now  introduce  a  state  feedback  variable,  v  (k),  according  to 

v{k)  =  Z-^  [u{k)  -  Frx{k)]  or  (11) 


u{k)  =  Frx{k)  +  Zrv{k).  (12) 

Substituting  Equation  12  back  into  9  and  10  the  new  system  in  terms  of  the  state  feedback  variable 
v{k)  is 

x{k  +  1)  =  [A  +  BFr]  x{k)  +  BZrv{k)  (13) 

y{k)  =  [C  +  DFr]  x{k)  +  DZrv{k)  (14) 
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Realizing  this  in  terms  of  a  RCF 


yields  the  following  results 


(15) 

(16) 


y  =  NM~^u. 

Fr  is  then  chosen  so  that  A  +  BFr  has  stable  poles.  In  past  research  ([4]  ,  [6] )  these  poles  have 
been  chosen  to  be  at  the  origin,  yielding  Finite  Impulse  Response  (FIR)  behavior.  Here,  however, 
the  only  requirement  is  that  the  poles  are  stable,  i.e.  they  must  be  within  the  unit  circle. 


Solving  the  dual  problem 

M  N]^ 

where  the  poles  of  .A  +  LC  are  chosen  so  the  system  is  stable. 


(17) 


2.2.3  State  Prediction 

Before  the  reference  signal,  v,  can  be  input  to  the  system  from  the  optimizer,  it  is  necessary 
first  to  calculate  that  signal  from  known  information.  Expanding  upon  Equations  12,  13,  and  14, 
the  discrete  time  state  estimates  are  written  as 

x{k  +  1)  =  Ax{k)  +  Bu{k)  +  L  \y{k)  -  y{k)]  (18) 

but  the  estimated 

y{k)  =  Cx{k)  +  Du{k)  (19) 

and 

u{k)  =  FrX{k)  +  ZrV{k) 
yields 
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x{k  +  l)  =  [A  +  BFr  +  LC  +  LDFr]x{k)  +  [BZr  +  LDZr]v{k)-Ly{k).  (20) 
Equation  20  can  be  rewritten  for  any  time  step  in  the  future  as 

x{k  +  l  +  l)  =  F{l)x{k  +  l)  +  G{l)v{k  +  l)  +  H{l)y{k  +  l)  (21) 

I  =  0...P-1 


where 


(22) 


(23) 


(24) 


F(0')  —  A  BFf  +  LC  "1“  LDFj- 

F{m)  =  A  +  BFr  ,  m  =  1 .  ..I 

G{1)  =  BZr  +  LDZr 

H{0)  =  -L 

H{rn)  =  0  ,  m  =  1 . .  J 
The  reason  that  H{m)  =  0  is  because  there  are  no  future  measured  outputs,  only  predictions.  Now 
it  is  possible  to  arrange  equations  21, 22, 23  and  24  to  produce  a  vector  of  future  predicted  states.  At 

this  point  the  assumption  that  D  is  zero  for  all  physical  systems  is  imposed.  Rewriting  Equation  21 

x{k  + 1) 

:  I  =  Fx{k)  +  Gv{k)  +  G°°v°°{k)  Hy{k)  (25) 

x{k+p)  ^ 

where 

m 


F  = 


G  = 


F(1)F(0) 


n  m 

L  i=p-i 

G 

F{l)G 


n  m 

j=r 


G 


0 

G 


n  m 

j=r 


G 


0 

0 


G 


1 

1 _ 

G 

n  F{j) 

G 

ft  m 

j=p-i 

j=p-i 

G 
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G°°  = 


0 

G 


0 

0 


H  = 


F{r  +  1)G 

G 

r+1 

n  Hj) 

G 

F{r  +  2)G 

j=r+2 

’  r+1 

r+2 

n  m 

G 

n  m 

j=p-i 

j=p-l 

H{0) 

1 

0 

0 


G 


n  m 

j=p-i 


G 


Fil)H{0) 


m 


n  ^0')! 

\j=p-i 

From  Equation  12,  any  future  control  input  is  written  as 

rt(A:  1)  —  Frx(^k  -F  /)  “h  Zj~v{k  -F  t) 

I  =  0...q-l 

u(k  + 1)  can  now  be  written  in  vector  form  as 
u{k) 


}  =  Fux{k)  +  Gum  +  G^v°°{k)  +  Huy{k) 


where 


Fu  = 


u{k  +  q  —  l)  J 

Fr 

FrF{0) 


n  m 


(26) 


(27) 


(28) 
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Gu  = 


GOO  _ 

u 


Hu  = 


J=9-2 


From  the  cost  function,  Equation  3,  a  method  for  determining  Au  is  required.  Letting 

Au{k)  =  u{k)  —  u{k  —  1) 
the  following  relationships  are  formed. 

f  'I 


[  Au{k  +  q-l)  j 


- 

Zr 

0 

0 

FrG 

Zr 

0 

1 

2 

Fr 

n  m 

G  Fr 

n  F{j) 

G 

Zr 

j-r~2 

j=r-2 

1 

2 

Fr 

n  m 

G  Fr 

n  m 

G  ••• 

FrG 

j=zr-l 

j=r~l 

1 

2 

r 

Fr 

n  F{j) 

G  Fr 

n  Fij) 

G  Fr 

n  F{j) 

G 

- 

J=q-2 

J=q-2 

_j=q-2 

- 

0 

0 

0 

Zr 

0 

0 

FrG 

Zr 

0 

FrF{r  + 1)(? 

FrG 

0 

r+1 

r+2 

\ 

Fr 

n  F{j) 

G  Fr 

n  F{j) 

G  •••  Fr 

U  Fij) 

G 

- 

J=q-2 

_j=q-2 

J-q-2 

- 

■ 

0 

FrH{0) 

Fr 

n  FU) 

m 

(29) 


(30) 


}  =  FAx{k)  +  GAHk)  +  G^v°°{k)  +  HAy{k)  +  Iu{k  -  1)  (31) 
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where 


i^A  = 


Ga  = 


GT 


Fr 

FrF{0)-Fr 
FrF{l)F{0)  -  FrF{0) 


'  0  0  " 

Fr 

n  F{j)-  n  F{j) 

j=q-3  j=q-^ 

’  0  o' 

Fr 

n  F{j)-  n  F{j) 

j=q-2  j=q-S 

Zr 

0 

- 

FrG  —  Zr 

0 

FrF{l)G  -  FrG 

•  •  • 

0 

1  1 

Fr 

n  m-  n  m 

G  ••• 

j=r—2  j=rS 

'  1  1  ' 

Fr 

n  FU)-  n  FU) 

G 

Fj-G  —  Z-r 

j~r—l  j=r—2 

'  1  1  ' 

’ 

r  r 

Fr 

n  ^’0)-  n  m 

G 

Fr 

n  m-  n  m 

G 

J-q-S  j=q-^ 

11 

1 

00 

<03. 

11 

1 

’  1  1  ' 

r  r 

Fr 

n  FU)-  n  FU) 

G  •••  Fr 

n  F{j^-  n  FU) 

G 

_j=q-2  j=q-3 

J=q-2  j=q-3 

- 

0 

0 

Zr 

0 

FrG  —  Zr 

0 

FrF{r  +  1)G  -  FrG 

0 

r+l  r+1 

p  p 

Fr 

n  n  p(j) 

G  •••  Fr 

n  m-  n  m 

G 

j=q-2  j-q-3 

II 

1 

to 

Oc. 

II 

1 

05 

_ 

Ha  = 


0 

FrH{0) 

FrF(0)H{0)  -  FrH{0) 


n  FU)-  n  Fij) 

j=q-2  j=q-3 


HiO) 


(32) 
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-I 

0 

J  =  .  . 

0 

Note,  these  prediction  matrices  are  based  on  the  absolute  control  input,  u,  being  input  into  the  plant. 
Previous  works  [4]  and  [6] ,  have  modified  the  plant  to  accept  incremental  control  inputs,  Au. 


2.2.4  Quadratic  Programing  Problem 

To  solve  the  cost  function  of  Equation  3,  Matlab’s  quadratic  programing  algorithm,  qp.m, 
was  used.  The  routine  is  based  upon  the  following  relationship 

J{x)  =  min  I  + /^a; I  (33) 

subject  to  Ax  <  b 

The  equations  for  the  state  predictions,  incremental  control  inputs  and  system  constraints  must  be 
rewritten  to  agree  with  the  format  of  Equation  33.  A  new  vector  reference  signal,  v{k),  is  defined 
in  Equation  34  and  is  the  vector  over  which  the  cost  function  is  minimized. 

v{k)  =  [  v{k)'^  v{k  +  1)^  •  •  ■  v{k  +  r  —  1)^  (34) 


2.2.4.1  Quadratic  Programing  Cost  Function 


Expanding  the  MMO  cost  function  of  Equation  3  and  using  the  state  prediction  matrices 
J{k)  =  v{k)^Sv{k)  +  [  x{k)'^  y{k)'^  v°°{k)'^  s{k)'^  ]Tv{k) +ind  ( 


where  S,  T  and  ind  are  given  by 


S  =  G'^C'^RyCG  +  GIRuGa 


T  =  2 


CF 

CH 

CG°° 


RyCG-\-2 


ind  = 


{Fx{k)  +  G^v°°  +  Hy{k)Y  CFkyC  (Fx{k) 


RiGa 


+  G°°v°^  + 


Hy{k)) 


{PAxik)  +  +  ^Ay(fc))  C^RaC  (FAx{k)  +  G^v'^  +  HAy{k))  + 

s{kfRy  [s(A:)  -  2C  (Fx{k)  +  Hy{k)  + 
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Since  ind  is  independent  of  the  optimization  variable  v,  it  was  neglected.  Additionally,  C,  Ry,  and 


Ru  are  defined  as 

C 

=  diag  (C, . . 

■,C) 

(39) 

Ry 

=  diag(iij„. 

•  *  5^) 

(40) 

Ru 

=  diag  {Ru, . 

•  •  5  Ra) 

(41) 

The  far  future  reference  input,  v°°,  is  calculated  using  Lemma  5.1  in  [6] . 

(p-r) 

C {A  +  BFry-'^  BZrV°°  =  s°°  (42) 

The  purpose  of  calculating  v°°  is  discussed  in  the  Pole  Placement  for  System  Stability  and  Steady 
State  Error  section. 


2.2.4.2  Quadratic  Programing  Constraints 

Ideally  there  would  be  unlimited  control  power  and  the  system  would  be  able  to  handle  un¬ 
limited  outputs.  However,  all  physical  systems  have  limits  on  their  control  inputs  and  occasionally 
they  have  implied  limits  on  the  system  outputs.  For  aerospace  systems  control  inputs  are  limited 
to  their  maximum  and  minimum  position  limits  as  well  as  maximum  and  minimum  rates.  Various 
outputs  may  also  be  limited  based  upon  potential  damage  to  the  system  or  exceeding  human  factor 
constraints.  It  is  the  ability  of  MPC  strategy  to  handle  system  constraints  within  the  control  problem 
which  makes  it  appealing  for  highly  dynamic  systems.  System  constraints  are  expressed  across  the 
control  and  prediction  horizons  as 

"b  ^  A  i')  ^  ”b  0 

i  =  0 ...  —  1 

A  ^)  ^  'u(/iJ  -j-  i)  ^  '^TncLX^^  ”b  (45) 


i  =  —  l 


“b  i) 


^  xQz  “b  ^)  ^  ^maxi}^  ”b  ^) 
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i  =  l...p 

where  q  and  p  are  the  control  and  prediction  horizons.  System  outputs  can  be  determined  by  a 
linear  combination  of  states  using  Equation  10,  thereby  allowing  output  constraints  as  well  as  state 
constraints. 

Defining  the  system  constraints  to  match  Matlab’s  quadratic  programing  format  in  Equation 
33,  the  following  relationships  are  constructed  [4] 

LxAu{k)  <  l{k) 


where 


and 


Mfiu{k)  <  m{k) 
x{k  + 1) 

Nu  I  <  n{k) 

x{k  +  p) 


Lx  = 

Ml,  = 

K  = 


-h,  1 
h,,  . 

-h, ' 

Ik,p 


(44) 


(45) 


m{k)  =  [  1  (46 

'^max 

n{k)  =  \  . 

^max 

Using  the  prediction  matrices  of  Equations  26, 29,  and  32  the  following  relationships  are  defined. 


Dv{k)  <  Ec{k) 

LxGa 
D  =  Mi,Gu 

N^G 

■  -LxFa  -LxGt  -Lx&a  -La  /  0  0  ■ 

E  =  -M„F^  -M„G^  -M„Hu  0  0  /  0  (47) 

-NyF  -N^G°°  -N^H  0  0  0  / 
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c(fc)  =  [  x{k)'^  v°°{k)'^  yi.k)'^  u{k  —  lY  m{kY'  n{k)'^ 

2.3  Pole  Placement  for  System  StabUity  and  Steady  State  Error 

Recent  work  [4]  and  [6]  has  focused  on  placing  the  poles  of  observer  and  state  estimator 
controllers  .Equations  15  and  17,  at  the  origin  for  discrete  systems.  This  has  assured  mathematically 
that  the  system  would  reach  a  desired  setpoint  in  a  finite  amount  of  time.  However,  this  assumes 
that  all  states  are  modeled  and  the  plant  is  an  LTI  model.  Using  real  systems,  not  all  the  states 
are  modeled,  noise  sources  introduce  error  and  generally  systems  are  nonlinear  and  time  varying. 
Because  of  these  potential  problems,  controller  pole  placement  was  relaxed  from  the  placement  of 
discrete  poles  at  the  origin  (continuous  at  — oo)  to  within  the  unit  circle  (stable).  Because  of  this 
relaxation,  it  was  not  possible  to  guarantee  trajectory  tracking  of  a  LTI  model  within  a  finite  amount 
of  time.  However  it  was  mathematically  possible  to  reach  a  window  around  the  tracking  point  in  a 
finite  amount  of  time. 

From  reference  [6] ,  Lemma  5.1,  a  relationship  between  the  far  future  reference  input,  v°°,  and 
the  far  future  setpoint  trajectory,  s°°,  is  given  by 

2k 

BFry~^  BZrV°°  =  s°°.  (48) 

j=i 

This  is  the  result  that  was  used  in  Equation  42  to  solve  for  v°°.  Equation  48  is  based  upon  the  poles 
of  A  +  BFr  and  A  +  LC  being  zero.  However  in  this  research  the  discrete  poles  were  simply  made 
stable.  Therefore  it  was  necessary  to  constmct  a  relationship  between  controller  pole  placement, 
prediction/control  horizons,  and  setpoint  error.  Starting  with  Equation  5.36  of  reference  [6]  (shown 

in  Equation  49),  it  is  possible  to  determine  the  error  of  the  stable  system. 

’x{k  +  l)]  _  \  A  +  BFr  -LC 

xlk  +  l)\  ~  [  0  A  +  LC\[xlk)y 

^  r  A  +  BFr  -LC 
2^  0  A  +  LC 

j=l  L 

x{k  +  l)  =  Future  Estimated  States 
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x{k  +  l)  =  x{k  +  l)  —  x{k  +  l) 


(50) 


For  a  system  where  the  controller  poles  are  at  the  origin  and  letting  I  >  2k  [6]  yields 

n  2k 

=  0. 


A  +  BFr  -LC 
0  A  +  LC 


(51) 


Knowing  the  far  future  setpoint  trajectory,  s°°,  the  far  future  reference  input  is  calculated  using 


Equation  48.  For  a  simply  stable  system 

A  +  BFr  -LC 
0  A  +  LC 


as 


(52) 


oo. 


Since  it  is  impractical  to  let  Z  — >  oo,  a  heuristic  approach  to  determining  the  error  using  a  first  order 
SISO  system  was  derived  and  then  applied  to  the  MIMO  system. 

First  recall  the  following  definitions. 

p  =  prediction  horizon 

q  =  control  horizon 

r  =  optimization  horizon  (53) 

Ts  =  Discrete  lime  Step 

p  =  max(p,q) 

Also  let  the  dynamic  setpoint  trajectory  be  fixed  to  at  the  end  of  the  optimization  horizon,  r. 
Letting 

I  =  p-r  (54) 

provides  I  time  steps  for  the  error  to  reach  a  finite  value.  Let  the  error  be  defined  as 

err  =  1  —  W  (55) 

where 


w  =  User  defined  >i\5ndow  of  Steady  State  Error. 


(56) 
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For  a  first  order  continuous  system  of  the  form 


His)  = 


Pc 


S+Pc 


subject  to  a  step  input  the  error  is  defined  as 

err  =  6^°*- 

Setting  the  time  t  =  l*Ta  and  solving  for  the  continuous  pole 

In  (1  —  w) 


Pc  = 


in 


Now  recalling  the  relationship  between  continuous  and  discrete  pole  locations 

In  (2) 


Pc  = 


(57) 


(58) 


(59) 


(60) 


and  combining  with  Equation  59,  a  relationship  between  pole  placement,  prediction/control  hori¬ 
zons,  and  setpoint  error  is  found 


In  (1  —  ty)  _  In  (z) 
_ 


ITs 


In(l-w) 

e  I 


(61) 


z*  =  (1~'“^) 

which  is  analogous  to  the  MIMO  Equation  52.  As  an  example  z  —  0.82  given  w  —  98%  and  I  —  20. 
For  the  MIMO  system  Equation  61  was  used  as  a  rough  guide  to  place  the  poles  of  A  4-  BFr  and 
A  -f-  LC  such  that 

Fn  =  Initial  Pole  Location  or  Pole  closest  to  Unit  Circle  for  A-\-BFr 

Fr,  =  (62) 


Li  =  Initial  Pole  Location  or  Pole  closest  to  Unit  Circle  for  A  -I-  LC 

Li  =  (1-10^)".  (63) 

All  other  controller  poles  were  then  placed  closer  to  the  origin  than  their  respective  starting  locations. 
The  placement  of  these  poles  was  accomplished  using  Ackerman’s  Formula. 


18 


2.4  Tiining  Parameters 


Having  established  the  mathematical  development  of  the  MPC  strategy  it  is  beneficial  to  list 
the  various  tuning  parameters  in  Table  1. 


Table  1 .  Tuning  Parameters 


Ihning  Parameter 

Description 

Ra 

Control  Input  Weighting  Matrix 

Ry 

Tracking  Error  Weighting  Matrix 

A  -\-  BFr  Pole  locations 

State  Gains  of  Inner  Feedback  Loop 

A  +  LC  Pole  locations 

State  Estimator  of  Inner  Feedback  Loop 

Ts 

Discrete  lime  Step 

P,Q,r 

Prediction,  Control  and  Optimation  Horizons 

l^m^n 

Control  Rate  and  Limit  Constraints  and  System  Constraints 

s 

Commanded  Trajectory 

Currently  specific  relationships  between  system  tracking  and  required  control  power  are  not 
developed  in  terms  of  these  tuning  parameters.  Relationships  are  generally  developed  using  a  heuris¬ 
tic  approach  and  numerous  simulations. 
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Chapter  3  -  Plant  Models,  Trajectories  and  MFC  State  Space 

Formulation 

3.1  F-16  Model 

3.1.1  Linear  Model 

Generation  of  linear  models  was  performed  using  various  Matlab  function  calls  and  script 
files  listed  in  Appendix  A.  Several  of  these  script  files  and  function  calls  were  based  upon  the 
FORTRAN  code  developed  by  Dr.  B.L.  Stevens  [12]  and  transcribed  into  Matlab  format  by  Dr. 
Brad  Liebst  and  several  of  his  students  [8]  .  A  new  front  end  script  file,  Modgen.m  in  Appendix 
A,  was  written  to  automate  the  generation  of  linear  models.  Modgen.m  requires  a  user  defined 
input  matrix  given  in  the  script  file  casel.m,  whose  rows  correspond  to  different  flight  conditions 
at  which  linear  models  are  to  be  generated.  The  four  elements  in  each  row  correspond  to  Mach 
number,  bank  angle  in  degrees,  rate  of  climb  in  jtjsec  and  altitude  in  jt.  Modgen.m  then  calls  the 
Matlab  function  trimmerm  [3] ,  which  calculates  the  control  inputs  necessary  to  make  the  state 
derivatives  zero.  Next,  Modgen.m  calls  jacob.m  [2]  ,  which  calculates  the  linearized  state  space 
model  based  upon  the  trimmed  flight  condition  from  trimmerm.  Finally  Modgen.m  outputs  the 
state  space  matrices  and  equilibrium  values  into  a  file  unique  to  the  specified  flight  condition.  Both 
trimmerm  and  jacob.m  are  designed  to  use  any  nonlinear  aircraft  model  obeying  the  equation 

x  =  f{x,u).  (64) 

3.1.2  F-16  Nonlinear  Model 

The  nonlinear  Matlab  function  subfl6a.m  [13] ,  provides  the  x  vector  specified  in  Equation 
64.  The  code  utilizes  extensive  look  up  tables  based  upon  wind-tunnel  data  developed  by  NASA- 
Langley[12, 124] .  The  nonlinear  model  currently  has  time  invariant  physical  properties  which  are 
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listed  in  Table  2.  The  actuator  physical  limits  and  rates  are  listed  in  Table  3.  The  nonlinear  model, 
subfl6a.m,  has  limitations  which  are  important  to  remember  when  performing  simulations.  These 
limitations,  listed  in  Table  4,  are  based  upon  the  bounds  of  the  wind  tunnel  data  which  was  collected 
by  NASA-Langley.  Control  input  rate  and  position  constraints,  vectors  I  and  m,  are  taken  from  the 
information  in  Table  3. 


Table  2.  F-16  Physical  Properties 


Property 

Value 

Weight 

20, 500  lbs 

Moment  of  Inertia,  Jxx 

9, 496  Slug  *  fi^ 

Moment  of  Inertia,  Jyy 

55, 814  Slug  *  fi‘‘ 

Moment  of  Inertia,  Jzz 

63, 100  Slug  *  ft^ 

Moment  of  Inertia,  Jxz 

982  Slug  *  ft^ 

Span 

30.0  ft 

Area 

mJP 

Mean  Aerodyanmic  Center 

11.32  ft 

Table  3.  Control  Actuator  Data 


Actuator 

Deflection  Limit 

Rate  Limit 

Time  Constant 

Elevator 

±25.0° 

60° /sec 

0.0495  sec  lag 

Ailerons 

±21.5° 

80° /sec 

0.0495  sec  lag 

Rudder 

±30.0° 

120°/sec 

0.0495  sec  lag 

Throttle  No  Afterburner 

0...  0.7699 

0.0  sec  lag 

Throttle,  With  Afterburner 

0.77...  1.0 

100^ 

_ S£C _ 

0.0  sec  lag 

Table  4.  Nonlinear  F-16  Model  Limitations 


Parameter 

Minimum  Limit 

Maximum  Limit 

Throttle  Position 

0.0 

1.0 

Altitude,  h 

Oft 

50,000  ft 

Mach  Number 

0.0 

~  0.6 

Angle  of  Attack,  a 

-10° 

45° 

Angle  of  Side  Slip,  P 

0 

o 

CO 

1 

O 

O 

CO 

1 

The  throttle  rate  limit,  Aut,  was  set  to  an  arbitrary  large  value,  100^.  This  was  done  to 
simulate  a  near  instantaneous  throttle  input  response.  However,  the  engine  dynamics  are  modeled 
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with  an  appropriate  delay,  so  that  there  is  not  an  instantaneous  change  in  thrust.  Note  that  the  throttle 
input,  ut,  is  unitless  as  it  represents  a  percentage  of  engine  thrust. 

3.1.3  Plant  Block  Diagrams 

In  order  to  simplify  the  controller  block  diagram  and  assist  in  the  modular  design,  the  plant 
model  was  grouped  into  a  single  SiMULiNK  subsystem.  A  specific  plant,  linear  or  nonlinear,  with 
the  correct  number  of  inputs  and  outputs  could  easily  be  placed  into  the  complete  MPC  controller. 

3.1.3. 1  F-16  Linear  Block  Diagram 

Figure  1  shows  the  linear  model  configuration.  It  is  a  standard  state  space  configuration  except 
that  in  the  SiMULiNK  Discrete  Plant  block  C'  =  JandD  =  0.D  was  set  to  zero  because  physical 
systems  tend  to  have  an  inherent  delay  in  any  input  to  output.  Also  the  generation  of  the  linear 
models  by  trimmerm  calculates  D  =  0  for  states  which  do  not  involve  acceleration.  Setting  C  =  I 
provides  the  designer  the  flexibility  of  calculating  the  outputs  through  a  linear  combination  of  the 
states,  X,  by  using  Simulink  matrix  gain  blocks,  K.  These  outputs  are  then  multiplexed  in  the 
order  in  which  the  MPC  controller  is  designed. 
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Figure  1.  F-16  Linear  Model  Block  Diagram 

3.1.3.2  F-16  Nonlinear  Block  Diagram 

Figure  2  shows  the  SiMULiNK  block  diagram  used  for  the  nonlinear  model. 


23 


Figure  2.  Nonlinear  F-16  Model  Block  Diagram 

The  input  to  this  block  diagram  is  the  output  of  the  controller,  u,  which  is  fed  into  a  zero  order 
hold.  The  zero  order  hold  ensures  that  the  continuous  integration  performed  by  the  nonlinear  plant 
is  fed  a  discrete  input  for  each  time  step,  T®.  Since  u  is  a  perturbation  signal  and  subfl6a.m  operates 
on  total  inputs,  the  equilibrium  input  value,  uq,  is  added  to  u  before  being  input  into  the  Matlab 
function  block.  The  heart  of  this  diagram  is  the  Matlab  function  block  where  the  nonlinear  state 
derivatives  are  calculated  using  subf  16a.m.  The  state  derivatives  are  then  integrated,  in  the  integra¬ 
tor  block  i ,  using  a  Matlab  proprietary  5*^  order  Runge  Kutta  equation  solver.  The  total  states  are 
then  fed  back  into  the  Matlab  function  block.  The  initial  conditions,  in  the  integrator  block,  are 
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set  to  the  equilibrium  state  values,  Xeq^a  >  calculated  by  trimmer.m.  The  four  graphical  output  blocks 
xd  long,  xd  lat,  x  long  and  ;c  lat  are  used  to  display  the  lateral  and  longitudinal  states  during  sim¬ 
ulations.  The  specific  lateral  and  longitudinal  states  are  selected  using  the  matrix  gain  blocks,  K, 
preceding  these  graphical  blocks.  Note  that  the  longitudinal  state  outputs  are  perturbation  outputs 
since  Xeq^n  was  subtracted  in  the  Sum3  block.  Figure  2  also  includes  a  state  feedback  block,  Kfeed. 
This  block  was  included  in  the  design  to  allow  for  continuous  stabilizing  feedback  but  was  set  to 
zero  in  this  thesis.  As  in  the  linear  block  diagram.  Figure  1,  the  specific  outputs  were  calculated  by 
a  linear  combination  of  the  states  determined  by  the  Simulink  matrix  gain  blocks  Phi,  alt.  Veloc¬ 
ity,  Beta  and  Thetadot.  These  outputs  were  then  multiplexed  to  provide  the  correct  output  vector, 
y,  to  the  MFC  controller.  Finally  the  unit  delay  block  preceding  the  output  provides  the  controller 
with  the  current  output  vector  while  a  new  vector  is  being  calculated  by  the  continuous  integration. 

3.1.4  Initial  Conditions 

The  initial  conditions,  Thble  5,  for  all  the  simulations  were  chosen  to  simulate  realistic  flight 
conditions.  A  clean  F-16  has  an  optimum  climb  airspeed  of  ^  400  knots.  Since  the  nonlinear  model 
is  limited  to  a  Mach  number  of  «  0.6,  a  low  starting  altitude  was  chosen  to  achieve  an  airspeed  close 
to  400  knots. 


Table  5.  Initial  Flight  Conditions 


Parameter 

Value 

Mach 

0.6 

Altitude,  h 

100  ft 

\felocity 

390.6  knots  669.7  ftfsec 

Bank  Angle,  (j) 

0° 

Rate  of  Climb,  h 

0  ft/ sec 

C.G.  Location 

30%mac 

The  control  input  equilibrium  values  are  given  in  Equation  65  for  a  30%  c.g.  location. 

uo  =  [  0.262  -1.544°  0°  0°  (65) 
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where 


T 

W  =  [  Ut  Ugl  Uail  '^rud  ] 

Using  Modgen.m,  linear  models  were  generated  using  the  initial  conditions  in  Tkble  5  and  varying 
h.  The  corresponding  eigenvalues  are  given  in  Table  6. 


Table  6.  Eigenvalues  at  Mach  =  0.6,  Alt  =  100ft,  Bank  Angle  =  0  deg,  C.G.  =  0.30 


Longitudinal* 

Rate  of  Climb 

h  =  100/t/sec 

h  =  200/t/sec 

h  =  300/t/sec 

Short  Period 

-1.59  ±  1.99i 

-1.59±1.99i 

-1.59±1.99i 

-1.59±1.99i 

Phugoid 

-0.013  ±  0.055i 

-0.0091  ±  0.053Z 

-0.0021  ±  0.047? 

Engine 

-1.0 

-1.0 

-1.0 

-5.0 

Lateral* 

Dutch  Roll 

-0.54  ±  4.12i 

-0.54  ±  4.12i 

-0.55±4.1U 

-0.55  ±  4.10? 

Roll 

-4.96 

-4.97 

-4.98 

Spiral 

-0.0033 

+0.0039 

*  Eigenvalues  are  in  continuous  time  with  units  rad/ sec 

Note  that  the  eigenvalues  vary  slightly  as  h  varies.  However  at  h  =  300/i/sec,  there  is  a 
sudden  change  in  the  engine  dynamics.  This  represents  a  transition  from  military  power  to  after¬ 
burner.  This  discrete  jump  occurs  in  the  nonlinear  model,  subfl6a.m,  when  the  throttle  input  is  0.77 
or  greater.  To  avoid  problems  associated  with  this  discrete  jump  two  additional  constraints  were 
placed  upon  the  simulations.  First,  the  throttle  was  limited  to  =  0.765.  Second,  the  altitude 
trajectories  were  designed  for  h  <  290  ftfsec  based  upon  the  results  of  Table  7. 

Table  7.  Throttle  Position  for  Different  h  "V^ues 


h  =  290/t/sec 

'^throtle  —  0.764 

h  =  ^00  ft/ sec 

'^hrotle  ~  0.774 

C.G.  =  30% 


3.1.5  Actuator  Constraints 

As  discussed,  the  throttle  actuator  was  limited  to  =  0.765.  The  mdder  was  also  limited 
due  to  the  initial  flight  condition.  In  actual  flight,  the  rudder  is  seldom  used  beyond  a  degree  or 
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two  at  flight  conditions  other  than  landing  configurations.  Because  this  is  an  up  and  away  flight 
condition  (i.e.  no  flaps,  no  gear,  much  greater  than  Ig  stall  velocity),  the  rudder  was  limited  to 
±0.5°.  Table  8  shows  the  actual  constraints  which  were  used  in  simulations. 


Table  8.  MFC  Actuator  Constraints 


Actuator 

Deflection  Limit 

Rate  Limit 

lime  Constant 

Elevator 

±25.0° 

60° /sec 

0.0495  sec  lag 

Ailerons 

±21.5° 

Rudder 

±0.5° 

0.0495  sec  lag 

Throttle  No  Afterburner 

0...  0.765 

■Slim 

0.0  sec  lag 

Since  the  controller  is  designed  to  operate  on  perturbation  inputs,  it  is  necessary  to  redefine 
the  constraints  in  terms  of  the  physical  limits  and  equilibrium  values. 

I^max  —  '^max  '^0 

f^min  ~  '^min  ^0  (66) 

Since  the  equilibrium  condition  implies  ii  =  0,  there  was  no  modification  of  the  controller  input 
rate  limits,  Imax  or  ^min- 

3.2  Dynamic  Hrajectories 

The  MFC  controller  setup  that  was  used  in  this  thesis  assumes  that  there  was  a  future  control 
input,  v°°,  which  was  found  by  a  constant  future  setpoint  trajectory  s°°.  In  past  research  [4]  and 
[6] ,  this  future  setpoint  trajectory  was  simply  the  step  input  that  was  commanded.  However,  for  a 
dynamic  trajectory  a  method  of  determining  s°°  is  required.  The  method  chosen  was  to  let  s°°  = 
constant  after  the  optimization  horizon,  r,  has  passed. 

s°°{l)  =  s{k  +  r)  (67) 

r  ±  1  <  I  <  p 

Figure  3  displays  a  graphical  based  upon 

r  —  10 
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Trajectory 


P 


30 


(68) 


q  =  30. 


k  +  I  Time  Steps 


Figure  3.  Optimization  and  Far  Future  Trajectory 


The  dynamic  trajectories  are  generated  at  each  time  step  using  the  Matlab  function  file 
trajpnt.m.  Trajpnt.m  is  called  during  the  optimization  at  each  time  step.  It  is  envisioned  that  instead 
of  a  prescribed  trajectory  which  is  known  at  the  start  of  the  simulation,  pilot  modeled  inputs  could 
be  substituted  for  trajpnt.m. 


3.2.1  Bank  Angle 

Two  different  trajectories  were  used  for  the  desired  dynamic  bank  angle  trajectories.  The  first 


was  a  simple  half  sine  wave,  shown  in  Figure  4,  according  to 


t  tgg 

f  ^  fss 


(69) 
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(j)^x  =  Maximum  Desired  Bank  Angle. 

Since  roll  rate,  is  strongly  coupled  with  aileron  input  [9,  Chapter  5]  and  it  is  desirable  to  generate 
smooth  control  inputs,  a  second  trajectory  where  the  starting  and  ending  derivatives  are  zero  was 
produced  using  a  cosine  wave,  as  shown  in  Equation  70  and  Figure  5: 


Time  In  Seconds  Time  in  Seconds 

Figure  4.  Bank  Angle  Trajectory  -  Half  Sine  Figure  5.  Bank  Angle  Trajectory  -  Full  Co- 
Wave,  tss  =  lOsec  sine  Wave,  tss  =  lOsec 
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Tim©  in  Seconds 


Figure  6.  Altitude  Trajectory  -  Half  Cosine  Wave,  tss  =  20sec 


3.3  MPC  State  Space  Formulation 

The  following  three  sections  describe  in  detail  the  block  diagrams,  controller  pole  placement 
technique  and  the  adaptive  constraint  techniques  utilized  in  the  simulations. 

3.3.1  Block  Diagrams 

The  entire  controller  incorporating  the  plant  is  shown  in  the  Simulink  block  diagram  of  Figure 
7: 
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Figure  7.  Model  Predictive  Control  -  Overall  Simulink  Block  Diagram 

The  block  diagram  is  divided  into  five  main  areas.  The  first  is  the  actual  plant.  Shown  in  the 
loop  is  the  Simulink  icon  Non  Linear,  which  contains  the  nonlinear  F-16  model.  As  discussed, 
any  linear  or  nonlinear  plant  can  be  substituted  there  provided  the  individual  elements  of  the  plant 
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(72) 


output  vector  correspond  to  the  desired  trajectory.  As  configured,  the  plant  output  vector  is 

y{k)^[cf>  h  V  (3  qf 
where  units  used  are  degrees,  feet  and  seconds.  The  output  vector  in  Equation  72  can  be  changed. 
If  it  is  changed,  the  corresponding  ‘desired  trajectory’  vector,  s,  must  also  be  changed  accordingly 
in  the  Matlab  code  trajpnt.m. 

Continuing  around  the  loop  from  the  plant,  the  state  estimator  is  reached.  The  estimator  design 
is  based  upon  the  results  of  Chapter  2.  The  estimator  operates  on  the  plant  outputs,  y,  and  the  plant 
inputs,  u.  Its  state  space  representation  is  given  by 

G 

Note  the  gains  for  A+BFr  and  A+LC  were  chosen  based  upon  a  linearized  model  at  the  operating 
point  described  in  Table  5.  Figure  8  shows  the  block  diagram  used  to  realize  Equation  73.  The  three 
SiMULiNK  input  blocks  in_l,  in_2  and  in_3  correspond  to  a  vector  of  constants,  the  control  input 
vector  [  u  y  ]  and  a  user  definable  gain  scheduling  parameter.  The  two  Matlab  function  blocks, 
A  Matrix  and  B  Matrix,  call  the  Matlab  functions  Amatrix.m  and  Bmatrix.m.  These  function  calls 
are  used  to  select  the  LTI  state  space  matrices  A,  B,  C,  Fr  and  L.  The  output  of  the  state  estimator 
is  the  vector  x. 
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Figure  8.  Simulink  State  Estimator  Block  Diagram 

The  controller  was  designed  with  the  flexibility  to  perform  on-line  gain  scheduling  based  upon 
a  user  defined  input  parameter.  Examples  of  these  parameters  could  include  y,y,u  or  others.  When¬ 
ever  a  linearized  state  space  matrix  or  prediction  matrix  is  required  by  any  portion  of  the  MFC  con¬ 
troller,  a  Matlab  function  call  is  made  to  matget.m.  Currently  matget.m  returns  the  LTl  state  space 
matrices  and  prediction  matrices  which  were  calculated  based  upon  the  initial  conditions  of  Table  5. 

The  next  point  in  the  loop  is  the  gain  block,  Fr,  for  the  estimated  states,  x.  The  Matlab  func¬ 
tion  block  in  Figure  9  calls  the  routine  Frwhich.m,  which  determines  the  gain  Fr  through  matget.m. 
Then  the  estimated  states,  x,  are  multiplied  by  Fr  and  are  output  to  the  summing  block  preceding 
the  plant  input.  The  three  input  blocks,  in_l,  in_2  and  in_3,  correspond  to  a  vector  of  constants, 
the  state  estimates,  x,  and  a  gain  scheduling  parameter. 
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Mux 


Figure  9.  State  Feedback  Gains  -  Simulink  Block  Diagram 

The  outputs  of  the  plant,  state  estimator,  the  previous  control  input,  u{k~-l),  and  user  defined 
control  parameter  are  then  fed  into  the  Optimizer  block.  The  zero  order  hold  is  placed  at  the  input 
of  the  optimizer  to  ensure  that  the  optimizer  works  on  a  discrete  basis.  The  output  of  the  optimizer 
is  the  reference  signal  v.  This  is  fed  into  the  summing  block  and  a  perturbation  control  signal  is  fed 
into  the  plant.  The  code,  optimize.m  in  Appendix  A,  calls  matget.m  to  determine  the  ITI  state  space 
and  prediction  matrices.  Also,  adaptive  constraint  techniques  are  employed  within  optimize.m,  if 
selected. 

The  fifth  section  of  Figure  7,  located  at  the  upper  left  hand  comer,  is  comprised  of  the  various 
constants  used  throughout  the  MFC  controller.  These  constants  are  multiplexed  into  a  vector  and 
then  distributed  to  the  various  MFC  controller  sections. 
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3.3.2  Pole  Placement  Ibchnique 

As  presented  in  Chapter  2,  controller  pole  placement  was  not  restricted  to  the  origin  in  this 
thesis.  Using  Ackerman’s  Formula  [5,  page  497]  a  pole  placement  routine  was  written,  gainsS.m, 
listed  in  Appendix  A.  The  routine  requires  starting  locations,  and  Li,  for  poles  closest  to  the  unit 
circle  of  A  +  BF^  and  A  +  LC  and  the  discrete  time  step,  Tg.  The  routine  then  generates  a  vector 
of  pole  locations  in  the  continuous  domain  where  the  bounds  are  given  by  Equation  74.  The  vector 
is  the  length  of  the  number  of  controller  states,  k. 


|Trnm|  — 
|Pmaa;|  ~ 


ln(T;j 


In  {PerFn) 


(74) 


Per  =  User  Defined  Percentage 

Once  the  gain  matrices,  Fr  and  L,  have  been  generated,  the  discrete  pole  locations  are  recomputed. 
If  the  pole  closet  to  the  unit  circle  is  greater  than  Fr,  or  Li,  the  corresponding  starting  value  is 
multiplied  by  the  percentage  R  and  the  appropriate  gain  matrix  is  recalculated.  For  this  thesis  Per 
and  R  equalled  95%  and  98%. 


3.3.3  Control  Input  Constraint  Weighting  Matrices 

In  Chapter  2  the  control  input  constraint  weighting  matrices  were  defined  as  shown  in  Equation 


75. 


Lx 


-k,, 

4.  J 


-4.9 

—Ik,p 

Ik,p 


(75) 


However,  for  the  simulations  Equation  76  lists  the  actual  values  assigned  to  these  weighting  matri- 


Lx 


-4,«  4, (9-1) 

.  4,e  4, (9-1) . 
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ces. 


(76) 


Mp  = 


Nu  =  [],  Null  Matrix 

By  choosing  Ni,  to  be  a  null  matrix  system  states  are  unconstrained.  The  structure  of  Lx  and 
places  constraints  on  the  control  inputs  at  the  time  step,  k,  but  during  the  optimization  problem  future 
control  inputs  of  A:  +  i  for  /  =  1 ...  9  are  unconstrained.  The  choice  of  these  constraint  matrices  was 
designed  to  provide  an  aggressive  control  scheme  for  a  highly  maneuverable  aircraft. 


3.3.4  Adaptive  Constraint  'Kchniques 

Two  adaptive  constraint  techniques  were  derived  to  assist  the  optimization  routine  in  finding 
the  global  minimum  at  each  time  step.  Initial  simulations  involving  the  LTI  model  showed  a  rela¬ 
tionship  between  roll  rate,  and  aileron  input.  This  agreed  with  conventional  thinking  [9,  75-77] 
and  for  a  simplified  model  the  relationship  between  roll  rate,  and  perturbation  aileron  input,  ASa, 
is  given  by  [9, 154,  Equation  5.6] 

^  ^  Z^ASa  (77) 

L/p 

Using  a  linear  simulation,  bank  angle  was  commanded  to  follow  a  sine  wave  input  as  described  in 
Equation  69.  The  resulting  bank  angle  output,  was  differentiated  and  plotted  against  aileron 
control  inputs.  A  linear  fit  was  applied  to  the  data  and  Equation  78  was  the  result. 

i  =  (78) 

Ua  =  —0.060 

In  the  code,  optimize.m,  0  was  approximated  by  a  first  order  finite  difference  method.  Also  an 
estimated  error  bank  rate  and  an  approximated  rate  of  climb  were  all  calculated  as  given  in  Equation 


79. 
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(79) 
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•  Sh  (A:  +  1)  —  Sfi  (k) 

=  - T, - 

The  adaptive  aileron  constraints  were  then  calculated  based  upon  three  separate  conditions  listed  in 
Table  9. 


Table  9.  Adaptive  Aileron  Selection  Criteria 


Case 

hp 

Equation  Number 

1 

0 

0 

80 

2 

0 

#0 

81 

3 

^0 

^0 

82 

Based  upon  the  conditions  of  the  three  cases,  new  maximum  and  minimum  actuator  constraint 
limits  are  defined.  In  Equations  80  through  82,  the  prime,  and  notation  indicates 

that  a  new  value  for  the  constraint  has  been  assigned.  The  maximum  and  minimum  values  rumaxa 
and  rriminay  correspond  to  the  original  actuator  maximum  and  minimum  constraint  values. 

m'n^_  = 


-0.06  i  0.25 

+  —rn. 


100 

-0.06  •,  0.25 


maxa 


+ 


100 


m„ 


(80) 


rur, 


•  0  5 

m'min,  =  -0-06^p -f 


(81) 


,  ^  •  -0.06  i  0.5 

=  -0-06(^p  + +  Yoo”^max„ 

"imin.  =  +  +  (82) 

The  0.5%  buffer  was  added  to  prevent  over  constraining  the  optimization  problem. 

The  elevator  adaptive  constraint  technique  was  derived  in  a  similar  manner.  The  altitude  output, 
h,  was  commanded  to  follow  a  cosine  trajectory  as  defined  in  Equation  71.  The  second  derivative 
of  h{t)  was  then  taken  and  plotted  against  the  elevator  input  signal  and  a  linear  fit  of  the  data  yielded 
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the  following 


_  ^alt 

“  “15  ■ 


(83) 


A  2%  to  5%  buffer  was  then  added  to  this  heuristic  equation  to  derive  the  following  adaptive  elevator 
constraints 

,  _  h  2.0 

^maXsi  -^QQ^Tnaxei 

,  _  h  2.0 

limine,  - 


(84) 


The  buffer  of  2%  to  5%  was  selected  for  two  specific  reasons.  First  there  was  no  consideration  of 
throttle  input  in  deriving  Equation  83.  This  introduces  error  into  the  assumption  that  h  is  dependent 
solely  on  elevator  deflection  since  classical  performance  techniques  have  shown  that  h  is  heavily 
dependent  upon  excess  thrust,  thus  implying  that  throttle  inputs  are  important.  The  second  reason  for 
the  buffer  was  because  there  was  no  consideration  of  required  elevator  during  bank  angle  changes. 
One  potential  fix  to  the  neglected  bank  angle  affect  is  to  derive  a  linear  relationship  between  bank 
angle  and  elevator  input  and  then  add  that  relationship  to  Equation  83. 


3.4  Operating  Envelope  and  F-16  Views 

In  Appendix  C  are  figures  of  the  F-16  operating  envelope  and  the  planform  and  longitudinal 
views  of  the  F-16  [1] .  The  planform  and  longitudinal  views  of  the  F-16  were  included  merely  as  a 
reference  for  the  reader.  The  operating  envelope  shows  the  specific  excess  power,  P^,  curves  for  a 
clean  F-16  at  military  power.  This  chart  was  included  to  show  that  the  initial  conditions  provide  a 
point  which  is  well  within  the  operating  envelope  of  the  F-16.  This  is  important  as  the  trajectories 
designed  for  this  thesis  attempted  to  simulate  realistic  conditions  showing  that  the  MFC  strategy  is 
viable  for  an  aerospace  system.  If  the  initial  conditions  were  on  the  edge  of  the  operating  envelope, 
the  realism  of  the  simulations  would  be  in  question.  Note  that  for  a  steady  level  turn  of  60°  (2  p’s) 
at  a  Mach  number  of  0.6,  Ps  is  almost  300  ft/ sec.  This  implies  that  for  constant  altitude  bank 
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angle  trajectories  where  the  bank  angle  does  not  exceed  60°,  the  throttle  should  never  be  required 
to  exceed  the  0.77  setting  corresponding  to  military  power. 
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Chapter  4  -  Results  of  Simulations 

Before  presenting  the  results  of  the  simulations  it  is  important  to  understand  three  basic  con¬ 
cepts.  First,  the  MFC  controller  consisting  of  the  MMO  cost  function,  state  prediction  matrices 
and  system  constraints  used  for  all  simulations,  was  designed  using  an  LTI F-16  model  at  the  initial 
conditions  presented  in  Chapter  3.  Results  which  are  described  as  linear  in  nature  imply  that  the 
plant  of  the  entire  system  was  the  same  LTI  F-16  model  used  for  the  construction  of  the  MFC  con¬ 
troller.  Simulation  results  which  are  described  as  nonlinear,  imply  that  the  plant  of  the  system  is  the 
nonlinear  F-16  model  shown  in  Chapter  3  using  the  Matlab  code  subfl6.m. 

4.1  Initial  Validation 

Before  attempting  to  use  a  new  controller  on  a  new  aircraft  model,  a  validation  of  the  MFC 
controller  design  and  dynamic  trajectory  tracking  was  required.  As  much  of  the  ground  work  nec¬ 
essary  for  this  thesis  was  completed  by  reference  [4] ,  the  LTI  model  used  in  that  thesis  of  the  High 
Angle  Research  "Sfehicle  (HARV)  was  used  to  validate  the  MFC  controller.  For  a  complete  descrip¬ 
tion  of  the  HARV  model,  operating  point,  scaled  inputs  and  scaled  outputs,  the  reader  is  referred  to 
reference  [4] .  The  state  space  discrete  realization  is  generated  by  a  Matlab  script  file,  harv.m. 
The  state  space  matrices  generated  are  shown  in  Appendix  C.  Equation  85  lists  the  order  of  states, 
inputs  and  outputs  with  input  descriptions  identified  in  Table  10. 
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Table  10.  Validation  Control  Input  Definitions 


Control  Input 

Description 

Stvs 

Symmetric  Thrust  \bctoring 

^AS 

Symmetric  Aileron 

^ss 

Symmetric  Stabilator 

^LES 

Symmetric  Leading  Edge  Flap 

^TES 

Symmetric  Trailing  Edge  Flap 

6t 

Throttle 

Equation  86  describes  the  trajectory  commanded  for  the  MFC  controller  validation. 


s  = 


0 


f  1.0(l-COs(jf;t)),  t 
\  0  t. 


tss  <  lOsec 
tss  >  lOsec 
ss  <  20sec 
>  20sec 


(86) 


-J 

As  seen  in  Figure  10  the  simulation  yielded  the  expected  sine  and  cosine  shapes  of  Equation  86. 
Figures  11  and  12  show  the  corresponding  control  inputs.  From  Reference  [4] ,  the  position  con¬ 


straint  limits  of  the  thrust  vectoring,  6tvs,  were  set  to  ±1  unit.  Figure  11  shows  that  Stvs  is  an 
active  constraint.  Figure  12  also  shows  an  active  rate  constraint  on  Leading  Edge  flaps,  SleSj  since 
its  rate  limit  was  set  to  ±0.5^^^  in  reference  [4]  .  These  active  constraints  serve  to  validate  the 
correct  functioning  of  the  MFC  controller  design. 


Figure  10.  Longitudinal  State  Outputs  of  ITI HARV  Model  using  Dynamic  Trajectories 
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Figure  11.  Control  Inputs  (1-3)  of  ITI HARV  Figure  12.  Control  Inputs  (4-6)  of  LTI HARV 
Model  using  Dynanndc  Trajectories  Model  using  Dynamic  Trajectories 

4.2  Linear  Model  Results 

Using  the  discrete  pole  placement  technique  described  in  Chapter  3,  a  sample  of  the  pole  lo¬ 
cations  using  Fn  =  0.75  is  giv.en  in  Equation  87.  Descriptions  of  simulations  list  the  starting  pole 
locations  and  Li. 

eig{A  +  BFr)  =  [  0.750  0.746  0.742  0.737  0.733  0.729  0.725  . . .  (87) 

0.721  0.717  0.713 

Equation  88  provides  the  discrete  time  step,  weighting  matrices  and  output  vector  description  which 

were  used  for  all  linear  simulations.  The  aircraft  model  trim  conditions  were  given  in  Chapter  3, 

Table  5,  along  with  the  control  input  equilibrium  values  in  Equation  65. 

Ts  =  0.1  sec 

0.1  0  0  0 
0  10  0 

0  0  10  0 
0  0  0  100 

■  100  0  0 

Ry  =  0  1000  0 

[  0  0  10000 

y{k)  =  [(j>  h  VY 


Tlm0  in  seconds 


42 


4.2.1  Bank  Angle  Response 


Seven  separate  cases,  listed  in  Table  11,  were  simulated  using  the  linear  model  to  characterize 


the  system  response  to  different  bank  angle  commanded  trajectories.  Cases  1  and  2  served  to  validate 
the  assumption  presented  in  Chapter  3  that  rudder  deflection  constraints  may  be  severely  limited 
without  sacrificing  tracking,  provided  a  landing  task  is  not  being  simulated.  Cases  2  and  3  further 
validate  that  simply  stable  pole  placement  provides  adequate  tracking.  Cases  4  through  7  were  used 
to  identify  the  optimum  bank  angle  trajectory  shape.  Also  the  advantages  and  disadvantages  of  the 
adaptive  aileron  constraint  technique  developed  in  Chapter  3  were  identified.  No  linear  simulations 
are  presented  here  using  the  adaptive  elevator  constraint  technique  as  it  did  not  affect  any  of  the 
linear  results. 


Table  11.  Linear  Model,  Bank  Angle  Cases 


Case 

A* 

B* 

C* 

D* 

1 

Off 

0.75 

±30° 

60°  Step  Input 

2 

Off 

0.75 

±0.5° 

60°  Step  Input 

3 

Off 

0.50 

±0.5° 

60°  Step  Input 

4 

Off 

0.75 

±0.5° 

60°  Sin  Input,  tss  —  lOsec 

5 

On 

0.75 

±0.5° 

60°  Sin  Input,  tss  —  lOsec 

6 

Off 

0.75 

±0.5° 

60°  Cos  Input,  tss  —  lOsec 

7 

On 

0.75 

±0.5° 

60"^  Cos  Input,  tss  =  lOsec 

*A  -  Aileron  Adaptive  Constraint  On  or  Off 

*B  -  j4  +  BFrand  A  +  LC  Pole  Starting  Location 

*C  -  Rudder  Stops 

*D  -  Trajectory 

Comparing  Figures  13  and  15,  there  is  no  noticeable  difference  in  rise  time  or  overshoot  for 
a  60°  step  input  despite  the  difference  in  rudder  limit  constraints.  While  the  variation  in  steady 
state  perturbation  outputs  V  and  h  in  Figures  14  and  16  are  negligible  compared  to  their  respective 
trim  conditions,  the  use  of  more  restrictive  rudder  constraints  did  increase  damping,  overshoot  and 
settling  time,  thus  improving  longitudinal  output  tracking.  A  comparison  of  Figure  15  and  the  Case 
3  corresponding  output  (see  Appendix  D)  showed  that  the  more  restrictive  pole  placement  of  Case 
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3  did  not  noticeably  improve  bank  angle  response.  It  did,  however,  slightly  improve  velocity  and 
altitude  response,  but  the  steady  state  tracking  response  in  Figure  16  was  more  than  acceptable  with 

AVss  =  -0.0005% 

Ahss  =  -0.01%. 


Time  (secs)  Time  (secs) 


Figure  13.  Case  1  -  Bank  Angle  Output  Figure  14.  Case  1  -  \felocity  and  Altitude  Out¬ 
puts 


Time  (secs)  Time  (secs) 


Figure  15.  Case  2  -  Bank  Angle  Output  Figure  16.  Case  2  -  \blocity  and  Altitude  Out¬ 
puts 

The  use  of  more  restrictive  rudder  constraints  also  showed  positive  benefits  for  the  control 
inputs.  Figure  20  shows  there  was  less  longitudinal  control  power  required  for  the  more  restrictive 
radder  constraints  as  compared  to  Figure  17.  There  is  a  much  more  dramatic  decrease  in  required 
lateral  control  power  as  seen  in  Figure  20.  As  seen  before  in  the  validation  simulation  results.  Figure 
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20  shows  the  aileron  rate  constraint  is  active.  Additionally  Figure  20  shows  the  rudder  position 
constraint  is  also  active.  The  aileron  is  experiencing  rate  constraints  according  to 

^'^max  ~  '^max  *  Tg 

^Umin  ~  '^hnin  *  Tg  (89) 

Where  Umaxau  =  S0°/sec  from  Table  8  of  Chapter  3  and  =  0.1  sec  yields  the  maximum  variation 
in  aileron  control  input,  Au,  between  time  steps  is  8°.  The  ratcheting  of  the  aileron  control  inputs 
is  due  to  the  unconstrained  future  control  input  u{k  + 1),  where  /  =  1 . . .  g,  as  discussed  in  Chapter 
3.  This  ratcheting  is  not  considered  to  be  a  problem  as  step  inputs  are  not  realistic  trajectories 
and  the  results  are  presented  here  to  demonstrate  the  capabilities  of  the  MFC  controller.  If  this 
ratcheting  were  to  create  problems  one  could  redefine  the  constraint  matrices,  Lx,  M^,  and  Nu  to 
allow  constraints  at  future  time  k  +  l.  A  limitation  of  the  ITI  model  is  shown  in  Figures  17  and  20. 
One  can  see  the  steady  state  values  of  elevator  and  throttle  are  both  zero.  An  actual  aircraft  would 
require  an  increase  in  both  (positive  deflection  of  throttle  and  negative  for  elevator)  to  maintain 
constant  altitude.  As  stated  in  the  Introduction,  this  research  concentrated  on  realistic  outputs.  As 
will  be  seen  in  the  nonlinear  simulations,  the  expected  increase  in  throttle  and  elevator  are  achieved. 
The  placement  of  the  MFC  controller  pole  locations  closer  to  the  origin  was  simulated  in  Case  3. 
The  results,  in  Appendix  D,  did  not  provide  any  significant  benefits  over  Case  2. 


45 


0.5 
0.4 
0.3 
0.2 

§  0.1 
o. 

I  0 
5-0.1 
-0.2 
-0.3 
-0.4 
-0.5 

0123456789  10  0123456789 

Time  (secs)  Time  (secs) 

Figure  17.  Case  1  -  Longitudinal  Controls  Figure  18.  Case  1  -  Lateral  Controls 
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Figure  19.  Case  2  -  Longitudinal  Controls  Figure  20.  Case  2  -  Lateral  Controls 

In  Cases  4  and  5  a  comparison  was  made  between  the  use  of  absolute  aileron  constraints,  Case 

4,  and  the  adaptive  aileron  constraint  technique  developed  in  Chapter  3,  Case  5.  The  results  are 
presented  in  Appendix  D  due  to  the  unrealistic  boundary  conditions.  The  sine  wave  bank  angle  tra¬ 
jectory  imposed  an  instantaneous  jump  in  roll  rate,  which  is  unrealistic  in  an  actual  aircraft.  The 
simulations  did  provide  informative  results  as  to  the  ability  of  the  adaptive  aileron  constraint  tech¬ 
nique  to  handle  these  discrete  boundary  conditions.  A  comparison  of  Cases  4  and  5  showed  that  the 
adaptive  constraint  technique  reduced  smooth  system  tracking.  This  is  attributed  to  the  suspected 
reduction  of  stability  caused  by  the  adaptive  aileron  constraint  technique.  Also  the  required  con¬ 
trol  power  was  increased.  Despite  these  consequences  of  the  adaptive  aileron  constraint  technique. 
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simulation  results  are  still  acceptable.  As  will  be  shown  later,  as  more  demanding  trajectories  are 
commanded,  the  benefits  of  the  adaptive  constraint  techniques  outweigh  these  consequences. 

In  Cases  6  and  7,  the  bank  angle  was  commanded  to  follow  a  cosine  trajectory.  The  use  of  a 
cosine  trajectory  imposed  zero  roll  rate  boundary  conditions,  which  provided  a  realistic  bank  angle 
trajectory.  Due  to  the  cosine  trajectory,  there  was  little  difference  between  the  results  of  Cases  6 
and  7  shown  in  Figures  21  through  28.  Figure  23  shows  a  much  improved  bank  angle  overshoot 
as  compared  to  Case  5.  Despite  slightly  increased  control  power  usage  due  to  the  adaptive  aileron 
technique,  the  results  of  Case  7  are  more  than  acceptable. 


Figure  21.  Case  6  -  Bank  Angle  Output 
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Figure  23.  Case  7  -  Bank  Angle  Output 


Figure  22.  Case  6  -  \blocity  and  Altitude  Out¬ 
puts 
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Figure  24.  Case  7  -  \blocity  and  Altitude  Out¬ 
puts 


Figure  25.  Case  6  -  Longitudinal  Controls  Figure  26.  Case  6  -  Lateral  Controls 


Figure  27.  Case  7  -  Longitudinal  Controls  Figure  28.  Case  7  -  Lateral  Controls 

When  the  adaptive  aileron  constraint  technique  is  employed,  it  is  the  term,  defined  in  Equa¬ 
tion  79  of  Chapter  3,  which  causes  the  oscillatory  behavior.  As  a  stability  analysis  was  not  per¬ 
formed  on  this  term,  these  oscillations  are  expected.  Additionally,  the  relative  contribution  of 
to  the  overall  adaptive  technique  can  cause  the  system  to  go  unstable.  For  this  thesis  a  trial  and 
error  approach  was  used  to  develop  the  appropriate  contribution,  trading  decreased  stability  for 
increased  tracking  performance.  Future  research  using  adaptive  constraints  should  address  these 
stability  issues  directly. 
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4.2.2  Altitude  Response 

Three  cases  listed  in  Table  12  were  simulated  to  characterize  commanded  altitude  response. 
Since  the  constrained  rudder,  ±0.5°,  showed  positive  results  in  the  bank  angle  responses  and  rud¬ 
der  is  primarily  a  lateral  control  input,  it  was  once  again  constrained  to  ±0.5°.  Additionally  there 
was  little  noticeable  improvement  in  placing  the  poles  of  the  controller  closer  to  the  origin  so  Iheir 
placement  remained  at  Fr^  =  0.75  and  Li  =  0.75. 

Table  12.  Linear  Model,  Altitude  Cases 


Case 

A* 

B* 

C* 

D* 

8 

Off 

0.75 

±0.5° 

500/f  Step  Input 

9 

Off 

0.75 

±0.5° 

500  ft  Cos  Input,  tss  =  20sec 

10 

On 

0.75 

±0.5° 

500/f  Cos  Input,  tss  =  20sec 

*A  -  Aileron  Adative  Constraint  On  or  Off 

*B  -  A  ±  BFrand  A  +  LC  Pole  Starting  Location 

*C  -  Rudder  Stops 

*D  -  Trajectory 

Figures  29  through  32  show  that  the  altitude  response  is  at  least  initially  unstable  when  com¬ 
manded  to  a  step  trajectory.  The  simulation  response  is  unrealistic  as  it  demands  an  infinite  rate  of 
climb  at  the  first  time  step.  Since  the  control  inputs  were  constrained,  this  implies  that  there  may 
not  be  a  feasible  solution  to  the  commanded  step  trajectory. 


49 


Tlm0  (secs)  Time  (secs) 


Figure  29.  Case  8  -  Altitude  Output  Figure  30.  Case  8  -  \felocity  and  Bank  Angle 

Outputs 


Time  (secs)  Time  (secs) 

Figure  31.  Case  8  -  Longitudinal  Controls  Figure  32.  Case  8  -  Lateral  Controls 

Given  sufficient  time,  the  system  should  in  theory  stabilize  and  converge  since  the  linear  theory 

states  that  the  cost  function  is  monotonically  decreasing,  however  the  response  shown  in  Figure  29 

would  be  objectionable  to  any  pilot.  The  violent  elevator  oscillations  would  also  be  objectionable 

as  they  would  decrease  the  life  expectancy  of  the  elevator  actuators.  The  ratcheting  of  the  elevator 

input  can  also  be  attributed  to  the  unconstrained  future  control  inputs  u{k  + 1)  as  discussed  earlier. 

Because  of  these  concerns  a  realistic  cosine  trajectory  was  commanded  in  Cases  9  and  10  according 

to  Equation  71  of  Chapter  3. 

Figures  33  through  36  show  the  simulation  results  using  a  commanded  cosine  trajectory  and  the 

affect  of  using  absolute  aileron  constraints  versus  adaptive  constraints.  In  both  cases  altitude  steady 
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state  tracking  error  is  less  than  4%.  However  the  big  improvement  of  using  adaptive  constraints  was 
the  significantly  reduced  steady  state  bank  angle  error,  27°  in  Figure  34  and  0.5°  in  Figure  36.  The 
corresponding  control  inputs  for  Cases  9  and  10  are  shown  in  Figures  37  through  40. 


Figure  33.  Case  9  -  Altitude  Output  Figure  34.  Case  9  -  \felocity  and  Bank  Angle 

Outputs 


Figure  35.  Case  10  -  Altitude  Output  Figure  36.  Case  10  -  \felocity  and  Bank  Angle 

Outputs 
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Figure  37.  Case  9  -  Longitudinal  Controls 
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Figure  39.  Case  10  -  Longitudinal  Controls 


Figure  38.  Case  9  -  Lateral  Controls 


Figure  40.  Case  10  -  Lateral  Controls 


4.2.3  Combined  Dynamic  Response 

The  final  two  cases  of  the  linear  simulations,  listed  in  Table  13,  show  the  ability  of  the  con¬ 
troller  to  handle  lateral  and  longitudinal  dynamics  simultaneously.  Since  the  linear  model  is  nearly 
decoupled  into  longitudinal  and  lateral  motion,  the  resulting  outputs  of  Case  11  were  basically  the 
results  of  Cases  5  and  10  superimposed.  Case  12  was  the  result  of  the  superposition  of  Cases  7  and 
10. 
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Table  13.  Linear  Model,Combined  Dynamic  Cases 


Case  II  A*  1  B*  I  C*  D* 


11  On  0.75  ±0.5°  60°  Sin  Input,  tss  =  lOsec  |  500ft  Sin  Input,  tss  =  20sec 


12  On  I  0.75  I  ±0.5°  |  60°  Cos  Input,  tss  =  Wsec  \  500ft  Cos  Input,  tgs  =  20sec 
*A  -  Aileron  Adative  Constraint  On  or  Off _ 


*B  -  A  ±  BFr^nd  A  ±  LC  Pole  Starting  Location 


*C  -  Rudder  Stops  _ _ 


*D  -  Bank  Angle  Trajectory  _ 


*E  -  Altitude  Trajectory  _  _ 
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Figure  45.  Combined  Case  11  -  Longitudinal  Figure  46.  Combined  Case  11  -  Lateral  Con- 
Controls  trols 


Figure  47.  Combined  Casel2  -  Longitudinal  Figure  48.  Combined  Case  12  -  Lateral  Con- 
Controls  trols 


4.3  Comparison  of  Linear  and  Nonlinear  Models 

Before  proceeding  with  the  nonlinear  F-16  model,  it  was  desired  to  first  ensure  that  the  non¬ 
linear  model.  Figure  2,  was  working  correctly  and  identify  any  potential  problem  areas  of  using  a 
linearized  model  to  develop  the  various  matrices  of  the  MFC  strategy.  To  achieve  these  objectives, 
the  linearized  F-16  model  of  Figure  1  was  placed  in  the  loop  of  the  MFC  controller.  Figure  7.  The 
system  was  commanded  according  to  equation  90,  where  the  output  vector  y(k)  is  defined  in  Equa¬ 
tion  91. 


s  =  [52000' 


(90) 
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y 


(91) 


=  h  V  ^  0] 

The  control  inputs  to  both  the  linear  and  nonlinear  systems  are  shown  in  Figures  49  and  50.  As  seen 
in  Figure  5 1  there  is  excellent  correlation  between  the  linear  and  nonlinear  model  for  bank  angle,  cf). 
Also  there  was  good  agreement  of  the  side  slip  angle,  /),  response  in  Figure  52  and  the  pitch  rate,  0, 
response  in  Figure  53.  However  the  two  areas  of  concern  are  with  altitude  and  velocity  responses 
shown  in  Figures  55  and  54.  These  results  caused  concern  as  the  state  prediction  would  be  predicting 
the  incorrect  future  values  of  velocity  and  altitude  which  are  fed  into  the  optimizer.  The  robustness 
of  the  MFC  strategy  would  be  questioned.  A  comparison  was  also  made  with  the  c.g.  location  at 
35%.  These  results  are  in  Appendix  D  and  show  worse  correlation  between  altitude  and  velocity. 


Time  (secs)  Time  (secs) 

Figure  49.  Linear  and  Nonlinear  Comparison  Figure  50.  Linear  and  Nonlinear  Comparison 
-  Longitudinal  Control  Inputs  -  Lateral  Control  Inputs 
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Bank  Angle  (deg) 


Time  (secs)  Time  (secs) 


Figure  5 1 .  Linear  and  Nonlinear  Comparison  Figure  52.  Linear  and  Nonlinear  Comparison 
-  Bank  Angle  -  Side  Slip  Angle 


Time  (secs) 


Figure  53.  Linear  and  Nonlinear  Comparison  -  Pitch  Rate 
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Time  (secs) 


Figure  54.  Linear  and  Nonlinear  Comparison  Figure  55.  Linear  and  Nonlinear  Comparison 
-  \felocity  -  Altitude 


4.4  Nonlinear  Model  Simulation  Results 

As  with  the  linear  simulations,  the  nonlinear  simulations  followed  a  similar  path  for  validating 
that  the  MFC  controller  works  on  both  lateral  and  longitudinal  trajectories  as  well  as  combined 
trajectories.  The  time  step,  weighting  matrices  and  output  vector  used  for  all  nonlinear  simulations 
are  given  in  Equation  92. 

Ts  =  0.1  sec 

■  0.1  0  0  O' 

0  le4  0  0 

^  ~  0  0  10 

0  0  0  100  _ 

■  le2  0  0  0  O' 

0  le4  0  0  0 

Ry  =  0  0  le4  0  0  (92) 

0  0  0  le2  0 

0  0  0  0  le4 

y  =  [ci>  h  V  p  ef 

The  added  outputs,  /3  and  6,  were  used  to  assist  in  system  tracking.  By  adding  these  two  outputs, 
two  additional  tuning  parameters  were  added  to  the  MFC  problem,  specifically  their  corresponding 
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elements  in  the  weighting  matrix  Ry.  The  indirect  benefits  of  these  states  were  seen  in  the  ability 
to  reduce  aileron  and  elevator  oscillations  without  adversely  affecting  system  tracking. 

4.4.1  Bank  Angle  Response 

Building  upon  the  linear  results,  5  different  nonlinear  bank  angle  response  cases  were  simu¬ 
lated.  These  cases,  listed  in  Table  14,  were  used  as  a  basis  of  comparing  the  linear  and  nonlinear 
responses.  In  general,  the  nonlinear  bank  angle  responses  closely  resembled  the  linear  responses  for 
the  same  set  of  conditions.  This  was  expected  as  the  comparison  between  the  linear  and  nonlinear 
models  showed  excellent  correlation  of  the  lateral  outputs. 


Table  14.  Nonlinear  Model,  Bank  Angle  Cases 


Case 

IIS 

B* 

C* 

D* 

Off 

■nBa 

ESQ 

60°  Step  Input 

Off 

0.75 

es&h 

60°  Sin  Input,  tgs  =  lOsec 

15 

On 

ESQ 

60°  Sin  Input,  tss  =  lOsec 

16 

Off 

ESQ 

60°  Cos  Input,  tss  =  lOsec 

17 

On 

0.75 

esb 

60°  Cos  Input,  tss  =  lOsec 

*A  -  Aileron  Adative  Constraint  On  or  Off 

*B  -  A  -f  BFrSnd  A  -I-  LC  Pole  Starting  Location 

*C  -  Rudder  Stops 

*D  -  Trajectory 

Case  13  demanded  a  60°  bank  step  trajectory.  Comparing  Figure  56  with  the  linear  case.  Figure 

15,  the  two  bank  angle  output  trajectories  are  very  similar.  Even  the  lateral  control  inputs.  Figures 

59  and  20,  show  similar  responses.  The  only  difference  is  the  nonlinear  case  took  approximately 

twice  the  time  to  reach  steady  state.  This  is  most  likely  due  to  aerodynamic  effects  not  accounted 

for  in  the  LTI  model.  An  interesting  result  of  the  nonlinear  simulation  was  the  steady  state  elevator 

and  throttle  inputs  shown  in  Figure  58.  As  discussed  in  the  linear  response  section,  the  steady  state 

longitudinal  control  inputs  were  zero.  Here  the  controller  is  simulating  a  real  world  response  in 

which  additional  elevator  and  throttle  inputs  are  required  to  generate  the  extra  lift  and  overcome 

the  additional  drag  generated  in  a  steady  level  turn.  While  the  controller  did  not  do  a  great  job  in 
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maintaining  zero  altitude  loss  shown  in  Figure  56,  a  40/t  drop  in  altitude  is  realistic  in  actual  flight 
for  a  sudden  60°  bank  angle  input.  Correct  modification  of  the  weighting  matrices  or  employing  an 
adaptive  weighting  matrix  scheme  should  reduce  this  altitude  deviation. 


Figure  56.  Case  13  -  Bank  Angle  and  Altitude  Figure  57.  Case  13  -\blocity,  Beta  and  Theta 
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Figure  58.  Case  13  -  Longitudinal  Control  In-  Figure  59.  Case  13  -  Lateral  Control  Inputs 
puts 

Simulation  results  for  Cases  14, 15  and  16  are  shown  in  Appendix  D.  For  Case  14  the  controller 
smoothly  followed  the  sine  wave  bank  trajectory,  but  lost  40  ft  in  altitude  during  the  dynamic 
portion.  Upon  reaching  zero  bank  angle,  the  altitude  output  also  returned  to  zero.  It  is  suspected, 
as  discussed  earlier,  that  the  inability  of  the  controller  to  accurately  predict  velocity  and  altitude 
far  future  outputs  is  the  cause  of  the  altitude  deviation.  Both  Cases  15  and  16  experienced  similar 


altitude  changes  before  returning  to  zero  deviation.  Once  the  output  trajectory  had  returned  to  zero 
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steady  state,  the  nonlinear  model  outputs  and  control  outputs  also  smoothly  returned  to  zero.  In 
Case  15,  as  discovered  in  the  linear  Case  5,  when  adaptive  aileron  constraints  were  utilized,  there 
was  not  as  smooth  of  tracking  for  the  same  commanded  trajectory  as  Case  14.  Aileron  inputs  were 
also  oscillatory  for  Case  15  but  followed  the  general  trend  as  Case  14. 

In  Case  16  the  bank  angle  response  was  similar  to  the  results  of  the  linear  Case  6.  Upon  reach¬ 
ing  the  zero  steady  state  value  at  10  sec  the  overshoot  was  negligible  especially  when  compared 
to  the  sin  wave  trajectory  of  Case  14.  This  reduced  overshoot  was  attributed  to  the  boundary  con¬ 
ditions  of  a  cosine  wave  as  discussed  in  the  linear  bank  angle  results  section.  Simulation  results 
for  Case  17  are  shown  in  Figures  60  through  63.  As  expected  from  earlier  results  the  trajectory  is 
not  as  smooth  since  adaptive  aileron  constraints  are  being  used.  However  there  is  still  acceptable 
trajectory  following  and  the  altitude  deviation  is  slightly  less  than  Cases  14  through  16. 


Figure  60.  Case  17  -  Bank  Angle  and  Altitude 


Figure  61.  Case  17  -  \blocity.  Beta  and  Theta 
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Figure  62.  Case  17  -  Longitudinal  Control  In-  Figure  63.  Case  17  -  Lateral  Control  Inputs 
puts 


4.4.2  Altitude  Response 

A  series  of  altitude  commanded  trajectories  were  simulated  as  listed  in  Table  15.  Cases  18 
through  20  were  run  to  compare  the  nonlinear  response  to  the  linear  results.  Based  upon  the  success 
of  the  Cases  18  -  20,  a  more  demanding  altitude  change  was  simulated  as  seen  in  Cases  21  through 
23.  Throughout  the  simulations  restrictive  radder  constraints  have  had  either  neutral  or  slightly  pos¬ 
itive  benefits  in  tracking  and  reduced  control  power.  Because  of  the  more  demanding  longitudinal 
maneuvers,  the  rudder  constraints  were  further  limited  to  ±0.1°.  Also  the  placement  of  the  state 
estimator  poles  with  Li  =  0.70,  appears  to  have  assisted  in  the  altitude  tracking.  However  quantifi¬ 
able  results  using  Li  =  0.70  were  not  identified. 
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Table  15.  Nonlinear  Model,  Altitude  Cases 


Case 

A* 

B* 

C* 

D* 

18 

Off 

0.75  and  0.75 

500/f  Step  Input 

19 

Off 

0.75  and  0.75 

20 

On 

0.75  and  0.75 

±0.5° 

21 

Off/Off 

0.75  and  0.70 

±0.1° 

9000ft  Cos  Input,  tss  —  60sec 

22 

On/Off 

0.75  and  0.70 

H- 
o 
i— i 
0 

9000/f  Cos  Input,  tss  =  60sec 

23 

On/On 

0.75  and  0.70 

±0.1° 

9000ft  Cos  Input,  tss  =  60sec 

*A  -  Aileron  Adative  Constraint  On  or  Ofl 

*B  -  A  +  BEpand  A  +  LC  Pole  Starting  Location 

*C  -  Rudder  Stops 

*D  -  Trajectory 

Due  to  instability  seen  in  the  linear  case  of  a  step  altitude  commanded  trajectory  Figure  29,  it 
was  informative  to  see  the  response  of  the  nonlinear  model.  Figures  64  through  68  .  The  nonlinear 
model  exhibited  much  better  system  tracking  as  compared  to  Case  8.  In  the  nonlinear  case,  the 
target  altitude  of  500  ft  was  not  reached  until  ~  7  sec  as  compared  to  3  sec  in  the  linear  simulation. 
It  is  expected  that  this  ‘sluggishness’  in  the  nonlinear  model  actually  increases  the  likelihood  of 
reaching  a  steady  state  value.  This  ‘sluggishness’  is  due  to  a  lack  of  understanding  of  the  limits  of 
various  accelerations  of  the  nonlinear  model  and  was  highlighted  in  other  simulations  not  included 
here.  For  example,  a  1000  ft  cosine  altitude  trajectory  with  tss  =  6  sec  was  commanded  of  the  LTI 
model.  While  the  commanded  rate  of  climb  only  reached  ~  250  ft/ sec  meeting  the  h  and  throttle 
requirements  of  Chapter  3,  there  was  no  consideration  of  the  aircraft’s  acceleration  response.  As  a 
result  the  nonlinear  model  went  unstable.  The  elevator  ratcheting  seen  in  Figure  67  was  suspected 
to  be  a  result  of  both  an  unrealistic  step  trajectory  and  the  unconstrained  control  inputs  u{k  +  1)  as 
discussed  in  the  linear  results  section. 
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Figure  64.  Case  18  -  Altitude  Output  Figure  65.  Case  18  -  \61ocity  and  Bank  Angle 

Outputs 
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Figure  66.  Case  18  -  Beta  and  Theta  Outputs 
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Figure  67.  Case  18  -  Longitudinal  Controls 


Figure  68.  Case  18  -  Lateral  Controls 


Case  19  (see  Appendix  D)  and  Case  20,  Figures  69  through  73,  have  near  identical  model 
output  and  control  input  results.  The  exception  is  that  the  use  of  adaptive  aileron  constraints  in  Case 
20,  proved  to  reduce  bank  angle  output  deviation  by  more  than  50%  over  the  first  20  sec.  In  both 
cases  though  steady  state  values  of  bank  angle  were  both  ~  0.4°. 

As  discussed  in  the  linear  bank  angle  results  section,  there  was  no  stability  analysis  made  on  the 
adaptive  aileron  constraint  technique.  This  proved  to  be  important  because  the  use  of  the  technique 
in  Equation  80  of  Chapter  3  proved  to  cause  the  output  of  Case  20  to  initially  become  unstable.  To 
regain  lateral  stability  the  buffer  term,  which  was  added  to  the  term  was  increased. 

This  provided  the  optimizer  more  latitude  in  selecting  a  stabilizing  aileron  input.  Equation  93  was 
the  original  form 


^maxa 

and  Equation  94  has  the  increased  buffer. 
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Figure  69.  Case  20  -  Altitude  Output  Figure  70.  Case  20  -  \felocity  and  Bank  Angle 

Outputs 
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Figure  71.  Case  20  -  Beta  and  Theta  Outputs 


Figure  72.  Case  20  -  Longitudinal  Controls  Figure  73.  Case  20  -  Lateral  Controls 

The  results  of  Cases  21  through  23  are  very  similar  to  each  other.  Plots  for  Cases  21  and  22  are 

found  in  Appendix  D  and  for  Case  23  Figures  74  through  78  are  seen  below.  The  use  of  the  adaptive 

elevator  constraint  developed  in  Equation  84  of  Chapter  3  had  no  apparent  affect  on  the  results  of 

the  simulations.  However  the  use  of  the  adaptive  aileron  constraint  technique  was  beneficial  for 

ensuring  bank  angle  deviation  remained  near  zero.  In  Case  21  the  bank  angle  reached  a  maximum 

of  10°  and  a  steady  state  value  of  8°.  In  Cases  22  and  23  the  maximum  bank  angle  was  less  than  1° 

as  well  as  the  steady  state  value.  As  discussed  and  shown  previously  the  only  draw  back  to  using 

adaptive  aileron  constraints  was  an  increase  in  aileron  control  power  due  to  oscillations.  However 
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Outputs 


as  seen  in  Figure  78  the  oscillations  are  negligible  when  compared  to  the  absolute  aileron  control 
inputs  of  ±21.5°. 


Time  (secs) 


Time  (secs) 


Figure  74.  Case  23  -  Altitude  Output  Figure  75.  Case  23  -  \felocity  and  Bank  Angle 

Outputs 
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Figure  76.  Case  23  -  Beta  and  Theta  Outputs 
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Figure  77.  Case  23  -Longitudinal  Inputs 


Figure  78.  Case  23  -  Lateral  Inputs 


4.4.3  Combined 

The  simulation  cases  listed  in  Table  16  were  run  to  meet  the  of  the  initial  objective  of  this 
thesis.  That  was  to  show  the  MFC  strategy  can  handle  combined  lateral  and  longitudinal  dynamics 
using  a  nonlinear  model.  Cases  24  and  25  were  simulated  to  compare  the  affects  of  adaptive  aileron 
constraints  during  the  combined  lateral  and  longitudinal  maneuvers.  Case  25  was  the  nonlinear 
equivalent  to  the  linear  Case  12.  Cases  26  through  28  were  run  to  demonstrate  the  MFC’s  robustness 
in  handling  commanded  trajectories  which  are  realistic  of  actual  aerospace  systems. 


Tkble  16.  Linear  Model,Combined  Dynamic  Cases 


Case 

A* 

B* 

C* 

D* 

E* 

'BHI 

Off/Off 

0.75  and  0.75 

±0.5° 

60°  Cos  Input,  tss  =  lOsec 

bOOft  Cos  Input,  tss  =  20sec 

On/Off 

0.75  and  0.75 

±0.5° 

60°  Cos  Input,  tss  =  lOsec 

500/f  Cos  Input,  tss  =  20sec 

26 

0.75  and  0.70 

±0.1° 

60°  Cos  Input,  tss  =  20sec 

9000/i  Cos  Input,  tss  —  60sec 

27 

On/Off 

0.75  and  0.70 

±0.1° 

60°  Cos  Input,  tss  =  20sec 

28 

On/On 

0.75  and  0.70 

±0.1° 

60°  Cos  Input,  tss  =  20sec 

9000/f  Cos  Input,  tss  =  60sec 

*A  -  Aileron  Adative  Constraint  On  or  Ofl 

*B-  A  +  BFrSnd  A  +  LC  Foie  Starting  Location 

*C  -  Rudder  Stops 

*D  -  Bank  Angle  Trajectory 

*E  -  Altitude  Trajectory 
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The  results  of  Case  24  (Appendix  D)  are  similar  to  Case  25  results  shown  in  Figures  79  through 
83.  In  both  cases  altitude  response  lagged  while  the  bank  angle  was  near  its  maximum  value.  With 
correct  initial  weighting  or  even  an  adaptive  weighting  scheme,  this  problem  should  be  solved.  The 
only  dilference  between  the  two  cases  was  the  usual  aileron  input  oscillations  caused  by  the  adaptive 
aileron  constraints.  The  adaptive  aileron  constraints  did  not  provide  any  beneficial  results  in  Case 
25  and  small  increase  in  control  power  was  required  as  compared  to  Case  24. 


Figure  79.  Case  25  -  Altitude  Output 


Figure  80.  Case  25  -  Bank  Angle  Output 


Figure  81.  Case  25  -  \felocity.  Beta  and  Theta  Outputs 
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Figure  82.  Case  25  -  Longitudinal  Controls  Figure  83.  Case  25  -  Lateral  Controls 

However,  the  use  of  adaptive  aileron  constraints  was  beneficial  for  system  tracking  as  more 

demanding  trajectories  were  commanded.  In  Case  26  (results  are  in  Appendix  D)  bank  angle  track¬ 
ing  was  smoother  than  in  Figure  85  of  Case  27  and  the  control  input  oscillations  seen  in  Figure  88 
were  not  experienced  in  Case  26.  But,  the  steady  state  bank  angle  in  Case  26  was  8°  and  it  was  less 
than  1°  in  Case  27.  The  other  outputs  and  inputs  shown  in  Figures  84  through  88  were  nearly  iden¬ 
tical  between  the  two  cases.  While  the  altitude  goal  of  9000/f  was  not  achieved,  its  steady  state 
error  was  less  than  3.5%.  With  additional  work  this  error  should  be  reduced.  Case  28  had  similar 
results  to  Case  27  as  seen  in  Figures  89  through  93.  The  difference  between  the  two  cases  was  the 
use  of  elevator  adaptive  constraints.  Figures  87  and  92  show  the  adaptive  elevator  constraints  trun¬ 
cated  the  elevator  input  at  10  sec.  Despite  this  truncation  there  was  no  apparent  affect  on  lateral  or 
longitudinal  outputs.  The  truncation  problem  should  be  overcome  by  simply  readdressing  the  ele¬ 
vator  adaptive  constraint  technique  and  including  a  lateral  contribution  as  discussed  in  Chapter  3. 
The  elevator  constraint  technique  was  included  in  this  thesis  to  demonstrate  it  had  no  adverse  affect 
on  system  tracking. 
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Figure  84.  Case  27  ^  Altitude  Output 
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Figure  85.  Case  27  -  Bank  Angle  Output 


Figure  86.  Case  27  -  \felocity.  Beta  and  Theta  Outputs 


Figure  87.  Case  27  -  Longitudinal  Controls  Figure  88.  Case  27  -  Lateral  Controls 


Figure  89.  Case  28  -  Altitude  Output  Figure  90.  Case28  -  Bank  Angle  Outputs 


Figure  91.  Case  28  -  \felocity.  Beta  and  Theta  Outputs 


Figure  92.  Case  28  -  Longitudinal  Controls 


Figure  93.  Case  28  -  Lateral  Controls 


Another  unique  feature  of  Cases  26  through  28  was  that  the  final  location  of  the  controller 


poles  had  moved  significantly.  Equation  95  lists  the  starting  and  ending  locations  of  the  poles. 


eig{A  +  BFr)start 


eig{A  +  BFr)end 


eig(^A  -1-  LC^  start 


eig{A  +  LC)end 


[  0.750  0.746  0.742  0.737  0.733  0.729  . . . 

0.725  0.721  0.717  0.713  ] 

[  0.886  ±  0.0853*  0.858  0.795  ±  0.127*  0.795  ±  0.089*  . . 
0.718  ±  0.241*  0.618  ] 

[  0.700  0.696  0.692  0.688  0.684  0.680  . . . 

0.677  0.673  0.669  0.665  ] 

[  0.919  0.796  ±  0.039*  0.730  0.721  ±  0.065*  . . . 

0.706  0.617  0.527  0.520  ] 


(95) 
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Chapter  5  -  Conclusions  and  Future  Research 


5.1  Conclusions 

As  shown,  the  original  objectives  of  using  the  MPC  strategy  were  achieved.  Numerous  sim¬ 
ulations  demonstrated  the  ability  of  the  MPC  strategy  to  handle  dynamic  trajectories  while  using 
a  nonlinear  model.  Also,  the  placement  of  the  controller  poles  at  the  origin  to  guarantee  system 
tracking  was  relaxed  to  a  simply,  stable  requirement  with  a  heuristic  guarantee  of  a  finite  error  us¬ 
ing  a  linear  model.  Finally  two  adaptive  constraint  techniques  were  developed  and  simulated.  The 
aileron  adaptive  constraint  technique  successfully  demonstrated  the  ability  to  reduce  lateral  steady 
state  error  with  minimal  increases  in  required  control  power.  The  adaptive  elevator  technique  was 
not  analyzed  in  depth  and  showed  only  neutral  effects  on  tracking  and  control  power  usage. 

Throughout  the  research,  an  attempt  was  made  to  identify  the  influence  of  the  various  tuning 
parameters  listed  in  Chapter  2  on  system  tracking  and  required  control  power.  However,  generaliza¬ 
tions  regarding  the  affects  of  these  tuning  parameters  on  system  performance  were  not  developed 
because  a  detailed  sensitivity  analysis  was  not  performed.  Additionally  the  influence  of  the  tun¬ 
ing  parameters  on  the  system  performance  was  suspected  to  be  nonlinear.  This  nonlinear  behavior 
was  demonstrated  through  various  simulations  as  individual  changes  were  made  in  the  placement 
of  controller  poles,  weighting  matrix  elements,  time  step  selection  and  dynamic  trajectory  design. 

Several  simulations  were  performed  in  which  controller  poles  were  placed  either  closer  to  or 
further  away  from  the  origin.  However  there  appeared  to  be  no  consistency  in  improved  system 
tracking  or  reduced  control  power.  In  general,  there  was  an  inverse  relationship  between  elements  of 
the  weighting  matrices,  Ry  and  Ru,  provided  there  was  a  strong  coupling  of  the  corresponding  inputs 
and  outputs  (i.e.  aileron  deflection  and  roll  rate).  However  there  were  cases  where  this  relationship 
did  not  hold  true.  Selection  of  the  time  step,  Tg,  also  demonstrated  nonlinear  effects.  Reduction  of 
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the  time  step  from  0.1  sec  to  less  than  0.05  sec  generally  made  the  system  unstable.  Increasing  the 
time  step  also  produced  unstable  results.  Dynamic  trajectory  design  was  also  an  important  factor  in 
achieving  stable  simulations,  \brious  polynomial  bank  angle  and  altitude  trajectories  were  designed 
with  smooth  boundary  conditions,  but  none  were  as  successful  as  the  sine  and  cosine  trajectories. 
Thus  more  complex  trajectory  design  might  be  required  to  consist  of  a  sum  of  a  finite  number  of  sine 
and  cosine  functions.  Also  for  dynamic  trajectory  design,  it  is  important  to  consider  the  physical 
limits  of  the  aircraft.  For  example,  a  trajectory  was  designed  where  the  rate  of  climb,  h,  fell  within 
the  bounds  of  the  nonlinear  model.  However  the  derivative  of  h,  was  beyond  the  capability  of  the 
nonlinear  model  and  caused  the  simulation  of  the  entire  closed  loop  system  to  become  unstable. 

5.2  Future  Research/Projects 

With  the  eventual  goal  of  realizing  the  MFC  strategy  on  an  actual  aircraft,  five  potential  re¬ 
search  topics  not  covered  in  this  thesis  are  presented.  First,  as  was  seen  in  Chapter  4,  linear  and 
nonlinear  states  quickly  diverged  from  each  other  when  provided  with  the  same  control  inputs.  It  is 
expected  that  a  simple  nonlinear  model  which  incorporates  the  effects  of  lateral  motion  would  solve 
this  problem.  The  prediction  matrices  would  then  be  rewritten  to  incorporate  this  nonlinear  model. 
A  simpler  interim  solution  would  base  the  linear  prediction  matrices  upon  an  ITI  model  which  was 
calculated  at  some  small  bank  angle.  However,  the  problem  of  lateral  and  longitudinal  coupling 
would  still  exist.  Second,  the  design  philosophy  used  in  this  thesis  separated  the  commanded  tra¬ 
jectory  code  from  the  optimization  code.  Doing  this  allows  the  incorporation  of  a  pilot  model  to 
determine  the  commanded  trajectory.  One  concept  of  a  simplified  pilot  model  would  consist  of  three 
parameters.  First,  longitudinal  stick  inputs  would  correspond  to  pitch  angle  changes.  Second,  lat¬ 
eral  stick  inputs  would  correspond  to  roll  rate  changes.  And  last  throttle  inputs  would  correspond  to 
altitude  changes.  The  third  idea  consists  of  developing  a  gain  scheduling  technique.  Since  the  dy- 
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namics  of  a  nonlinear  model  are  constantly  changing,  the  use  of  gain  scheduling  could  significantly 
improve  system  tracking.  The  design  used  in  this  thesis  left  open  the  possibility  of  selecting  the  pre¬ 
diction  matrices  as  well  as  the  inner  loop  feedback  matrices.  The  selection  criteria  could  be  based 
upon  any  combination  of  outputs,  control  inputs  or  their  derivatives.  The  selection  would  take  place 
at  each  time  step  within  the  optimization  routine.  Fourth,  to  simulate  more  realistic  conditions,  the 
SimStar  hybrid  computer  is  available  here  at  ART  Using  this  asset,  the  cumbersome  requirement 
of  numerical  integration  would  be  performed  by  the  analog  portion  of  the  SimStar  computer.  And 
finally,  perhaps  the  most  difficult  of  these  topics  is  the  identification  of  the  effects  of  the  various 
tuning  parameters  on  system  performance. 
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APPENDIX  A  -  MatlabCode 


A.l  Modgen.m 

% 

%  in  Model  Generation  for  Thesis  using  Aircraft  Control 
%  and  Simulation  by  Stevens  and  Lewis 
%  Capt  Christopher  M.  Shearer 
%  1  Jul  97 

% 

clear 

global  ay  az 

xguess  =  [669  0.01  0  0.  0.01  0  0  0  0  0  0  100  16.98]; 
xguess  =  [669.796  0.0111544  0  -6.27372e-10 ... 

0.0111544  0-0-0  0  0  0... 

100  16.9845]; 


uguess  =  [0.261541  -1.54463  6.18088e-ll  -2.02108e-ll]; 
const(l)  =  0.0;  %  const(l)  =  gamma  (deg) 

%  const(2)  =  Sin(Gamma) 
const(3)  =  0.0;  %  const(3)  =  Roll  Rate  (deg/sec) 
const(4)  =  0.0;  %  const(4)  =  Pitch  Rate  (deg/sec) 
const(5)  =  0.0;  %  const(5)  =  Turn  Rate  (deg/sec) 
const(6)  =  0.0;  %  const(6)  =  Phi  (Bank  Angle)  (deg) 

%  const(7)  =  Cos(Phi) 

%  const(8)  =  Sin(Phi) 

const(9)  =  0.0;  %  const(9)  =  Theta  dot  (deg/sec) 

const(lO)  =  1;  %  const(lO)  =  Coordnated  Turn  (1  =  yes,  0  =  no) 

const(ll)  =  0.0;  %  const(ll)  =  Rate  of  Climb  (ft/sec) 

const(12)  =  0.0;  %  const(12)  =  Bleed  Rate  (ft/sec'^2) 

const(13)  =  5  %  const(13)  =  Orient  (1  =  Wings  Level,  Steady  Flight) 

%  (2  =  Constant  Climb  Angle,  Gamma) 

%  (3  =  Constant  Altitude,  Constant  Tim  Rate) 

%  (4  =  Constant  Pitch  Rate,  theta  dot) 

%  (5  =  Constant  Turn  Rate  and  Constant  Climb) 
cases  =  easel 
n  =  size(cases,l); 


76 


for  i  =  l:n 
mach  =  cases(i,l) 
alt  =  ca 
ses(i,2); 

TR  =  cases(i,3); 

RC  =  cases(i,4); 

fname  =  [’vm’,num2str(mach*10),’a’,num2str(alt/100),’tr',... 

num2str(TR*  10);rc’  ,num2str(RC/100),\out^]; 

fid  =  fopen(fname/w’); 

xguess(12)  =  alt; 

xguess(l)  =  add  (mach, alt); 

const(5)  =  TR; 

const(ll)  =  RC; 

%%%%%  Go  find  steady  state  values  of  states  %%%%% 

[x,u,fcost,lcost]  =  trimmer(xguess,uguess,const,fid); 

%%%%%  Go  compute  the  state  space  matrices  %%%%% 

[a,b,c,d]  =  jacobl(x,u); 

%%%%%  Get  the  X  dot  terms 
xd  =  subfl6(x,u); 

[amach,qbar]=adc(x(  1 )  ,x(  1 2)) ; 

%%%%%  Output  the  results 

dataout(x,u,const,xd,fcost,lcost,ay,az,qbar,amach,a,b,c,d,fid); 

%%%%%  Close  the  current  file 
fclose(fid); 

%%%%%  Set  the  current  steady  state  values  to  the  guess  for  the  next  time 
%  xguess  -  x; 

%  uguess  =  u; 
end 
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A.2  TtimmeKm 

function  [x,u,fcost,lcost]=trimmer(Xguess,Uguess, const, fid) 

% 

%  Program:  trimmer 
%  by:  Mech  628  Incredible  Group  1 
%  Ise,  Shearer,  Clark 

% 

%  This  program  numerically  calculates  the  equilibrium  state  and  control  vectors  of  an  F- 16  model  given 
%  certain  parameters.  Inputs  include  initial  guesses  for  the  equilibrium  state  and  input  vectors. 

%  If  the  routine  is  called  with  no  inputs  the  user  will  be  prompted  to  key  the  equilibrium  initial 
%  guesses  in  by  hand. 

%  The  user  will  be  prompted  to  pick  one  of  the  following  A/C  orientation  options 
%  and  provide  the  desired  altitude,  airspeed,  gamma,  turn  rate,  pitch  rate, etc.  : 

% 

%  1.  Wings  Level  (gamma  =  0) 

%  2.  Wings  Level  (gamma  <>  0) 

%  3.  Steady  Constant  Altitude  TUrh 
%  4.  Steady  Pull  Up 
%  5.  Steady  Tim  and  Constant  Climb 
% 

%  The  user  will  also  be  prompted  for  the  number  of  iterations  to  be  used  in  the  numerical 
%  minimization  search. 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  states:  controls: 

%  xl  =  Vt  x4  =  phi  x7  =  p  xlO  =  pn  ul  =  throttle 
%  x2  =  alpha  x5  =  theta  x8  =  q  xll  =  pe  u2  =  elevator 
%  x3  =  beta  x6  =  psi  x9  =  r  xl2  =  alt  u3  =  aileron 
%  xl3  =  pow  u4  =  mdder 
% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Script/Function  calls: 

%  getinput 
%  subfl6 
%  clfl6 
%  confl6 
%  fminsa 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
global  ay  az 
format  long; 
if(nargin=:=4) 
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x=:Xguess; 

u=Uguess; 

else 

x=zeros(l,13); 

u=zeros(l,4); 

end 

% 

%  const(l)  =  gamma  (deg) 

%  const(2)  =  Sin(Gamma) 

%  const(3)  =  Roll  Rate  (deg/sec) 

%  const(4)  =  Pitch  Rate  (deg/sec) 

%  const(5)  =  Turn  Rate  (deg/sec) 

%  const(6)  =  Phi  (Roll  Angle)  (deg) 

%  const(7)  =  Cos(Phi) 

%  const(8)  =  Sin(Phi) 

%  const(9)  =  Theta  dot  (deg/sec) 

%  const(lO)  =  Coordnated  Turn  (1  =  yes,  0  =  no) 

%  const(ll)  =  Rate  of  Climb 
%  const(12)  =  Bleed  Rate 

%  const(13)  =  Orient  (1  =  Wmgs  Level,  Steady  Flight) 
%  (2  =  Constant  Climb  Angle,  Gamma) 

%  (3  =  Constant  Altitude,  Constant  Turn  Rate) 

%  (4  =  Constant  Pitch  Rate,  theta  dot) 

%  (5  =  Constant  Turn  Rate  and  Constant  Climb) 

% 

rtod  =  180/pi; 
orient  =  const(13); 
ndof  =  6; 

const(l)  =  const(l)/rtod; 
const(2)  =  sin(const(l)); 
const(3)  =  const(3)/rtod; 
const(4)  =  const(4)/rtod; 
const(5)  =  const(5)/rtod; 
const(6)  =  const(6)/rtod; 
const(7)  =  cos(const(6)); 
const(8)  =  sin(const(6)); 
const(9)  =  const(9)/rtod; 
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clear  s 

if  (orient  =  3)  |  (orient  ==  5) 
s(l)=u(l); 
s(2)=u(2); 
s(3)=u(3); 
s(4)=u(4); 
s(5)=x(2); 
s(6)=x(4); 
s(7)=x(5); 
end 

if  (orient  ==  1)  |  (orient  ==  2)  |  (orient  ==  4) 
s(l)=u(l); 
s(2)=u(2); 
s(3)=x(2); 
end 

cont=  100.0; 
itertot  =:  0.0; 
iter=  1500; 
tol  =  0.5; 

%%%%%  Now  solve  for  the  trimmed  condition  as  long  as  the  change  in 
%%%%%  the  cost  function  is  greater  than  0.5  percent 
while  cont  >  tol 

options  =  [0  l.OE-9  l.OE-9 0000000000 iter]; 

[s, options, x,u,fcost,lcost]  =  fminsa(’clfl6’,s,options,[],x,u,const); 

fprintf(’Initial  Cost  Function:  %g\n\  fcost); 

fprintf(Tinal  Cost  Function:  %g\n’,  Icost); 

cont  =  abs((fcost-lcost)/fcost)  *  100 

if  Icost  <  l.Oe-5 

cont  =  tol*0.9; 

end 

itertot  =  itertot  +  iter 
end 
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A.3  Jacob.m 

function  [a,b,c,d]=jacob(Xequil,Uequil) 
%  [a,b,c,d]=jacob(xequil,uequil) 


%  Program :  Jacob 
%  by:  Mech  628  Incredible  Group  2 
%  Capts  Chapa  &  St.  Germain 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  This  program  numerically  calculates  the  linearized  A,  B,  C,  &  D  matrices  of  an  F- 16  model  given 
%  certain  parameters.  Inputs  include  the  equilibrium  state  vector  and  input  vector.  If  the  routine 
%  is  called  with  no  inputs  the  user  will  be  prompted  to  key  the  equilibrium  values  in  by  hand. 

%  Output  may  be  shown  coupled  longitudinal/lateral  or  separate.  If  convergence  is  not  reached  for 
%  any  value,  the  user  is  prompted  for  an  estimate. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  states:  controls: 

%  xl  =  Vt  x4  =  phi  x7  =  p  xlO  =  pn  ul  =  throttle 
%  x2  =  alpha  x5  =  theta  x8  =  q  xll  =  pe  u2  =  elevator 
%  x3  =  beta  x6  =  psi  x9  =  r  xl2  =  alt  u3  =  aileron 
%  xl3  =  pow  u4  =  rudder 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Order  of  the  measurements  described  by  y  =  Cx  +  Du  are: 

% 

%  Longintudinal  (Elevator  Only):  y  =  [  Az  q  alpha  theta  Vt  ]T 
%  Longintudinal  (Elevator  and  Throttle):  y  =  [  Az  q  alpha  theta  Vt  ]T 
%  Lateral:  y  =  [  Ay  p  r  beta  phi  ]T 

%  Coupled  Long  &  Lat:  y  =  [  Az  q  alpha  theta  Vt  Ay  p  r  beta  phi  ]T 
% 

%  Note:  angles  are  in  degrees,  angular  rates  are  in  deg/s, 

%  airspeed  in  ft/sec,  accelerations  (Az,  Ay)  are  non-dimensional  (i.e.,  in  g’s) 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Order  of  the  inputs  u  are: 

% 

%  Longintudinal  (Elevator  Only):  u  =  [  el  ]T 
%  Longintudinal  (Elevator  and  Throttle):  u  =  [  thtl  el]  T 
%  Lateral:  u  =  [  ail  rdr  ]T 
%  Coupled  Long  &  Lat:  u  =  [  thtl  el  ail  rdr  ]T 
% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  Script/Function  calls: 
%  getinput 
%  subfl6 


global  az  ay; 
if  nargin==2 
x=Xequil; 
u=Uequil; 
else 

dispC  ’) 

dispC Input  The  Equilibrium  State  And  Control  \fectors:’) 
dispC  ’) 
getinput 
end 

xe=x;ue=u; 

tol=0.0001; 

xde=subfl6(x,u); 

aze=az;aye=ay; 

%%%%%  A  matrix  %%%%% 
for i=l:13 
for  j=l:13 

x=xe;del=0.01;slopel=0;diff=.9; 

if  xe(i)==0 

del=0.5; 

else 

del=del*xe(i); 

end 

while  diff  >  tol 

x(i)=xe(i)+del; 

xd=subfl6(x,u); 

slope2=(xd(j)-xde(j))/(del); 

diff=abs(slopel -slope2); 

del=:del*.l; 

slope  l=slope2; 

ifdiff>le6 

fprintf(’No  convergence  for  perturbed  state  XI  with  state  X2  :  ’); 

Xl=i 

X2=j 

slopel=input(’Enter  estimate :  ’); 

diff=0; 

end 

end 

AA(j,i)=slopel; 

end 

end 
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%%%%%  B  matrix  %%%%% 

x=xe; 

for  i=l:4 

forj=l:13 

u=ue;del=0.01  ;slopel=0;diff=  1 ; 

if  ue(i)=:=0 

del=0.5; 

else 

del=del*ue(i); 

end 

while  diff  >  tol 

u(i)=ue(i)+del; 

xd=subfl6(x,u); 

slope2=(xd(j  )-xde(j  ))/(del) ; 

diff=abs(slope2-slopel); 

del=del*.l; 

slopel=slope2; 

ifdiff>le6 

fprintf(’No  convergence  for  perturbed  input  U  with  state  X  :  ’); 

U=i 

X=j 

slope l=input(’ Enter  estimate  :  ’); 

diff=0; 

end 

end 

BB(j,i)=slopel; 

end 

end 
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%%%%%  C  matrix  %%%%% 
u=ue; 

for  i=l:13  %  az 

x=xe;del=0.01;slopel=0;diff=l ; 

if  xe(i)==0 

del=0.5; 

else 

del=del*xe(i); 

end 

while  diff  >  tol 

x(i)=xe(i)+del; 

xd=subfl6(x,u); 

slope2=(az-aze)/(del) ; 

diff=abs(slope2-slopel); 

del=del*.l; 

slopel=slope2; 

ifdiff>le6 

fprintf(’No  convergence  for  az  with  state  X  : 
X=i 

slopel=input(’Enter  estimate : 

diif=0; 

end 

end 

CC(l,i)=slopel/(-32.2); 

end 

for  i=l:13  %  ay 

x=xe;del=0.01  ;slopel=0;diff=l ; 

if  xe(i)==0 

del=0.5; 

else 

del=del*xe(i); 

end 

while  diff  >  tol 
x(i)=xe(i)+del; 
xd=subfl6(x,u); 
slope2=(ay-  aye)/ (del) ; 
diff=abs(slope2-slopel ); 
del=del*.l; 
slope  l=slope2; 
ifdiff>le6 

fprintf(’No  convergence  for  ay  with  state  X  :  ’); 
X=:i 
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slopel=input(’Enter  estimate :  ’); 

diff=0; 

end 

end 

CC(6,i)=slopel/(32.2); 

end 

rtod  =  180/pi; 

CC(2,:)=[0  000000  rtod  0000  0];CC(3,:)=[0  rtod  0  0  0  0  0  0  0  0  0  0  0]; 
CC(4,:)=[0  000  rtod  0000000  0];CC(7,:)=[0  0  0  0  0  0  rtod  0  0  0  0  0  0]; 
CC(8,:)=[0  0000000  rtod  000  0];CC(9,:)=[0  0  rtod  0  0  0  0  0  0  0  0  0  0]; 
CC(10,:)=[0 00 rtod 00000000 0];CC(5, :)=[!.  00000000000  0]; 
%%%%%  D  matrix  %%%%% 
x=xe; 

for  i=l:4  %  az 

u=ue;del=0.01;slopel=0;diff=l; 

if  ue(i)==0 

del=0.5; 

else 

del=del*ue(i); 

end 

while  diff  >  tol 

u(i)=ue(i)+del; 

xd=subfl6(x,u); 

slope2=(az-aze)/(del); 

diff=abs(slope2-slopel); 

del=deP.l; 

slope  l=slope2; 

ifdiff>le6 

fprintf(’No  convergence  for  az  with  input  U  :  ’); 

U=i 

slopel=inputCEnter  estimate :  ’); 

diff=0; 

end 

end 

D(l,i)=slopel/(-32.2); 

end 
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for  i=l:4  %  ay 

u=ue;del=0.01  ;slopel=0;diff=l ; 

if  ue(i)==0 

del=0.5; 

else 

del=del*ue(i); 

end 

while  diff  >  tol 

u(i)=ue(i)+del; 

xd=subfl6(x,u); 

slope2=(ay-aye)/(del); 

diff=abs(slope2-'Slopel); 

del=del*.l; 

slope  l=slope2; 

ifdiff>le6 

fprintf(’No  convergence  for  ay  with  input  U  : 

U=i; 

slopel=input(’Enter  estimate :  ’); 

diff=0; 

end 

end 

D(6,i)=slopel/(32.2); 

end 

D(2,:)=[0  0  0  0];D(3,:)=[0  0  0  0];D(4,:)=[0  0  0  0];D(7,:)=[0  0  0  0]; 

D(8,:)=[0  0  0  0];D(9,:)=[0  0  0  0];D(10,:)=[0  0  0  0];D(5,:)=[0  0  0  0]; 

%%%%%  transform  A,B»C  to  book’s  format  %%%%% 

posit=[l  25  8  13  3  47  9]; 
fori=l:9 

A(l,i)=AA(l,posit(i));A(2,i)=AA(2,posit(i));A(3,i)=AA(5,posit(i)); 

A(4,i)=AA(8,posit(i));A(5,i)=AA(13,posit(i));A(6,i)=AA(3,posit(i)); 

A(7,i)=AA(4,posit(i));A(8,i)=AA(7,posit(i));A(9,i)=AA(9,posit(i)); 

end 

B(1,:)=BB(1,:);B(2,:)=BB(2,:);B(3,:)=BB(5,:);B(4,:)=BB(8,:);B(5,:)=BB(13,:); 

B(6,:)=BB(3,:);B(7,:)=BB(4,:);B(8,:)=BB(7,:);B(9,:)=BB(9,:); 

C(:,1)=CC(:,1);C(:,2)=CC(:,2);C(:,3)=CC(:,5);C(:,4)=CC(:,8);C(:,5)=CC(:,13); 

C(:,6)=CC(:,3);C(:,7)=CC(:,4);C(:,8)=CC(:,7);C(:,9)=CC(:,9); 

%%%%%  output  to  a,b,c,d  %%%%% 
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dispC  ’) 

fprintf(’  1)  Longitudinal  set  (az,q,alpha, theta)  With  elevator  only  for  input  \n’); 
fprintf(’2)  Longitudinal  set  (az,q, alpha, theta)  With  elevator  and  throttle  for  inputs  \n’); 
fprintf(’3)  Lateral  set  (ay,pr,r,beta,phi)  \n’); 
fprintf(’4)  Combined  Lateral  and  Longitudinal  set  \n\n’); 

%choice2=input(’Choose  your  output  format :  ’); 

%%%%%  set  choice2  to  1  for  MECH  628  final 
choice2  =  4 
if  choice2==l 

A1=1;A2=4;B1=2;B2=2;C1=1;C2=5;D1=1;D2=5;C3=1;C4=4;D3=2;D4=2; 
elseif  choice2==2 

A1=1;A2=5;B1=1;B2=2;C1=1;C2=5;D1=1;D2=5;C3=1;C4=5;D3=1;D4=2; 
elseif  choice2==3 

A1=6;A2=9;B1=3;B2=4;C1=6;C2=10;D1=6;D2=10;C3=6;C4=9;D3=3;D4=4; 

else 

A1=1;A2=9;B1=1;B2=4;C1=1;C2=10;D1=1;D2=10;C3=1;C4=9;D3=1;D4=4; 

end 

a=A(Al:A2,Al:A2); 

b=B(Al:A2,Bl:B2); 

c=C(Cl:C2,C3:C4); 

d=D(Dl:D2,D3:D4); 
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A.4  Runme.m 

%  Orginal  author  Derek  W  Ebdon 
%  Modified  Extensively  by  Capt  Chris  Shearer,  Jul  97 


%  Script  file  for  controller  setup,  ’’setup. m” 
%  Lateral  and  Longitudinal  motion 
%  CONTROL  SURFACES:  All  functioning 


%  The  initial  flight  condition  is  as  follows: 

%  Alt.:  100  ft 
%  Mach:  0.6 
%Vt:  670  ft/s 
%  AOA  :  0.57  deg 
%  PA  :  0.57  (pitch  angle) 

%  FPA  :  0  deg  (flight  path  angle) 
clear 

bankstep  =  0; 
altstep  =  1 ; 
if  bankstep  ==  1 
Frpole  =  0.75; 

Lpole  =  0.75; 

Ts  =  0.1;  %  Time  Step 
end 

if  altstep  ==  1 
Frpole  =  0.75; 

Lpole  =  0.75; 

Ts  =  0.1;%  Time  Step 

ifTs==  0.051 

Frpole  =  0.85 

Lpole  =  0.85 

end 

end 

Tstart  =  0.0;  %  Start  Time 
Tstop  =  30.0;  %  Stop  Time 
%  Plant  States,  Inputs  and  Outputs 
%  x(l)  =  air  speed,  VT  (ft/sec) 

%  x(2)  =  angle  of  attack,  alpha  (rad) 

%  x(3)  =  pitch  angle,  theta  (rad) 

%  x(4)  =  pitch  rate,  Q  (rad/sec) 

%  x(5)  =  engine  thrust  dynamics  lag  state,  pow 
%  x(6)  =  angle  of  sideslip,  beta  (rad) 

%  x(7)  =  roll  angle,  phi  (rad) 

%  x(8)  =  roll  rate,  P  (rad/sec) 

%  x(9)  =  yaw  rate,  R  (rad/sec) 

%  x(10)  =  altitude,  h  (feet) 


88 


%  u(l)  =  throttle  command  0.0  <  u(l)  <  1.0 
%  u(2)  =  elevator  command  in  degrees 
%  u(3)  =  aileron  command  in  degrees 
%  u(4)  =  rudder  command  in  degrees 
%  y(l)  =  Az,  Normal  Acceleration  (g’s) 

%  y(2)  =  pitch  rate,  Q  (deg/sec) 

%  y(3)  =  angle  of  attack,  alpha  (deg) 

%  y(4)  =  pitch  angle,  theta  (deg) 

%  y(5)  =  air  speed,  VT  (ft/sec) 

%  y(6)  =  Ay,  Lateral  Acceleration  (g’s) 

%  y(7)  =  roll  rate,  P  (deg/sec) 

%  y(8)  =  yaw  rate,  R  (deg/sec) 

%  y(9)  =  angle  of  sideslip,  beta  (deg) 

%  y(10)  =  roll  angle,  phi  (deg) 

%  y(ll)  =  altitude  (feet) 

% 

%  Get  the  trimed  ITI  state  space  matrices 
% 

[ascl,bscl,cscl,dreal,dscl,xeq,xeqall,ueq]  =  trim3; 
eig(ascl) 

along  =  ascl(  1:5, 1:5); 

eig(along) 

alat  =  ascl(6:9,6:9); 

eig(alat) 

% 

%  Establish  the  maximum  and  minimum  control  inputs  for  the  model 
%  Note  these  are  in  deg  and  deg/sec,  page  584  of  Aircraft  Control 
%  and  Simulation 
% 

uabsmax  =  [0.765  25.0  21.5  0.5]; 
uabsmin  =  [0  -25.0  -21.5  -0.5]; 
umaxrate  =  [le2  60.0  80.0  120.0]; 
uminrate  =  [-le2  -60.0  -80.0  -120.0]; 

% 

%  The  new  output  states 
% 

cscl  =  [cscl(10,:);  %  Bank  angle  in  deg 
cscl(ll,:);  %  Altitude  in  feet 
100000000  0;%  \felocity  (feet/sec) 

0  0  0  0  0  180/pi  0  0  0  0;  %  beta.  Side  Slip  angle  (deg) 

0  0  0  180/pi  0  0  0  0  0  0];  %  theta  dot,  pitch  rate  (deg/sec) 
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dscl  =  [dscl(10,:); 
dscl(ll,:); 

zeros(l,size(bscl,2)); 
zeros(l  ,size(bscl,2)); 
zeros(l  ,size(bscl,2))]; 

% 

%  Discritize  the  model 
% 

[A,B,C,D]  =  c2dm(ascl,bscl,cscl,dscl,Ts,’zoh’); 

% 

%  Check  the  rank  of  the  observability  and  controllability 
% 

if  rank(ctrb(A,B))""=size(A,l)  |  rank(obsv(A,C))~=size(A,l) 
dispC’The  rank  of  the  Controllability  Matrix  is  ’); 
rank(ctrb(A,B)) 

disp(’The  rank  of  the  Observability  Matrix  is  ’); 
rank(obsv(A,C)) 

disp(’The  size  of  the  A  matrix  is’); 
size(A) 
end 
% 

%Set  the  continous  time  matrices 
% 

at  =  ascl; 
bt  =  bscl; 
ct  =  cscl; 
dt  =  dscl; 

% 

%Deterniine  the  number  of  control  inputs,  outputs,  and  states 
% 

kappa  =  size(A,l);  %estimator  states 
xi  =  size(B,2);  %control  inputs 
eta  =  size(C,l);  %system  outputs 
% 

%  Establish  the  prediction,  optimization,  and  control  horizons 
% 

p  =  26;  %  State  Prediction  Horizon 
q  =  26;  %  Control  Horizon 
r  =  5;  %  Optimization  Horizon 
if  min([p  q])  -  2*kappa  <=  r 

disp(The  Optimization  Horizon  is  within  2  times  the  number  of  states  to  min(p,q)’); 

break 

end 
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% 

rho  =  max([p  q]); 

% 

%  Establish  Z  matrix  weights 
% 

Zr  =  eye(xi); 

Z1  =  eye(xi); 

% 

%  Establish  R  optimization  weights 


if  bankstep  ==  1 
Ry  =  eye(eta); 

Ry(l,l)  =  le2*Ry(l,l);  %  Bank  Angle  (deg) 

Ry(2,2)  =  le4*Ry(2,2);  %  Altitude  (Ft) 

Ry(3,3)  =  le4*Ry(3,3);  %  Pert  \blocity  (ft/sec) 

Ry(4,4)  =  le2*Ry(4,4);  %  Side  Slip  Angle  (deg) 

Ry(5,5)  =  le4*Ry(5,5);  %  Pitch  Rate  (deg/sec) 

Ru  =  eye(xi); 

Ru(l,l)  =  le-l*Ru(l,l); 

Ru(2,2)  =  le4*Ru(2,2); 

Ru(3,3)  =  lel*Ru(3,3); 

Ru(4,4)  =  le2*Ru(4,4); 
end 

if  altstep  ==  1 
Ry  =  eye(eta); 

Ry(l,l)  =  le2*Ry(l,l);  %  Bank  Angle  (deg) 

Ry(2,2)  =  le4*Ry(2,2);  %  Altitude  (Ft) 

Ry(3,3)  =  le4*Ry(3,3);  %  Pert  \blocity  (ft/sec) 

Ry(4,4)  =  le2*Ry(4,4);  %  Side  Slip  Angle  (deg) 

Ry(5,5)  =  le4*Ry(5,5);  %  Pitch  Rate  (deg/sec) 

Ru  =  eye(xi); 

Ru(l,l)=le-l*Ru(l,l); 

Ru(2,2)  =  le4*Ru(2,2); 

Ru(3,3)  =  lel*Ru(3,3); 

Ru(4,4)  =  le2*Ru(4,4); 
end 
% 

%  Establish  the  constraint  matrices:  LI,  Mm,  Nn 

%  The  zeros  are  thrown  in  so  that  we  do  not  want  to  constrain  the  future  control 
%  inputs  but  rather  only  the  next  step  i.e.  the  identity  matrix.  This 
%  is  done  because  we  will  solve  for  the  control  inputs  down  the  line 
%  but  we  will  only  implement  the  first  step  set  of  control  inputs 
% 
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LI  =  [  -eye(xi)  zeros(xi,xi*(q-l)) ; 
eye(xi)  zeros(xi,xi*(q-l))  ]; 

Mm  =  [  -eye(xi)  zeros(xi,xi*(q-l)) ; 
eye(xi)  zeros(xi,xi*(q-l))  ]; 

Nn  =  []; 

% 

%Nn  =  [  »eye(kappa)  zeros(kappa,kappa*(p-l)) ; 

%  eye(kappa)  zeros(kappa, kappa* (p- 1 ))] ; 

% 

%Establish  physical  constraint  limits:  1,  m,  n 
% 

1=  l*[(-uminrate*I^)  (umaxrate*Ts)]’; 
m  =  l*[-(uabsmin-ueq)  (uabsmax-ueq)]’; 
n=[]; 

% 

%  Go  calculate  the  gains  neccesary  to  make  the  system  stable 
% 

[Fr,L,K]=gains3(A,B,C,D,'R,Frpole,Lpole); 

% 

%  Set  up  the  controller  matrices  based  upon  the  gains 
% 

[Ax,Bx,Cx,Dx,Ay,By,Cy,Dy,FO,Fl,G,H]=contller(A,B,C,Fr,L,Zl,Zr); 

% 

%  Establish  the  gain  matirices  for  the  non  linear  simulation 
% 

Kfeed  =  [K(l,l)  K(l,2)  K(l,6)  K(l,7)  K(l,3)  0  K(l,8)  K(l,4)  K(l,9)  0  0  K(l,10)  K(l,5); 
K(2,l)  K(2,2)  K(2,6)  K(2,7)  K(2,3)  0  K(2,8)  K(2,4)  K(2,9)  0  0  K(2,10)  K(2,5); 

K(3,l)  K(3,2)  K(3,6)  K(3,7)  K(3.3)  0  K(3,8)  K(3,4)  K(3,9)  0  0  K(3,10)  K(3,5); 

K(4,l)  K(4,2)  K(4,6)  K(4,7)  K(4,3)  0  K(4,8)  K(4,4)  K(4,9)  0  0  K(4,10)  K(4,5)]; 

Kfeed  =  zeros(4,13); 

P=:  180/pi; 

% 

%  Non  linear  Matrix  Gains,  Used  for  picking  off  outputs 
% 

Krc  =  [0000000000010]; 

Kh  =  Krc; 

Kphi  =  [0  00P00000000  0]; 

Kvel  =  [100000000000  0]; 

Kbeta  =  [0  0P0  00  00  0  00  00]; 

Kthetd  =  [OOOOOOOPOOOOO]; 
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Klookl  =  [1  0  00  0  00  00  00  00;  %\fel 
0P0000000000  0;%  alpha  (deg) 
0000P0000000  0;% Pitch  (deg) 
0000000P0000  0;%Q  Pitch  Rate  (deg) 
0000000000010;%  altitude  (feet) 
0000000000001];%  engine  thrust 
Klook2  =  [OOPOOOOOOOOOO;%  beta  (deg) 
OOOPOOOOOOOOO;%phi  (deg) 
OOOOOOPOOOOOO;%E  Roll  Rate  (deg) 
OOOOOPOOOOOOO;%Yaw  (deg) 

0  0  0  0  0  0  0  0  P  0  0  0  0];  %  R,  yaw  rate  (deg) 

% 

%  Linear  Matrix  Gains,  Used  for  picking  off  outputs 


KvelL  =  [1000000000]; 

KbetaL  =[00000P0000]; 

KphiL  =  [000000P000]; 

KqdotL=[000P000000]; 

KhL  =  [0000000001]; 

% 

%  Set  up  the  prediction  matrices 
% 

[cF,cG,cGinf,cH,cFu,cGu,cGuinf,cHu,cFdel,cGdel,cGdelinf,cHdel,cIdel] 

predmat(kappa,xi,eta,p,q,r,rho,FO,Fl,G,H,Fr,Zr); 

% 

%  Set  up  the  contraints 
% 

[S,T,cD,cE]  =  qpconst(C,Ry,Ru,cF,cG,cGinf,cH,cFdel,cGdel,cHdel, ... 
cGdelinf,cFu,cGu,cHu,cGuinf,Ll,Mm,Nn,p,q,r,rho); 

% 
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%  Establish  the  initial  values 


% 

xhatO  =  0.0*ones(l, kappa); 
xO  =  zeros(l, kappa); 
uO  =  zeros(l,xi); 

% 

%  Set  the  number  of  x  dot  terms  coming  out  of  the  xd  routine 
%  for  the  F-16  Model 
% 

kapfl6  =  13; 

x0fl6  =  zeros(l,kapfl6); 

% 

%  Save  the  Matrices  to  be  used  later 
% 

save  matsl  A  B  C  D  Fr  L  Zr  T  S  cD  cE  1  m  n  ueq 

ailu  =  []; 

aill  =  []; 

save  ail  ailu  aill 

% 

%  Start  the  Simulink  diagram  for  the  non-linear  simulation 
% 

nlinsim 
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A.5  GainsS.m 

function  [Fr,L,K]=gains3(A,B,C,D3,zstgoal,zesgoal); 


%  function  [Fr,L,K]=gains3(A,B,C,D,Ts,zstates,zestim) 
% 

%  Function  designed  to  make  the  poles  of  the  discrete 
%  system  A  +  B*Fr  less  than  zstates  and  the  poles  of  the 
%  discrete  system  A  +  L*C  less  than  zestim 
% 

%  Fr,  Discrete  Gains 
%  K,  Continuous  Gains 
% 

%  Convert  the  system  to  a  continous  models 
% 

[Ac,Bc,Cc,Dc]  =  d2cm(A,B,C,D,Ts,’zoh’); 


%  Get  the  number  of  states 
kappa  =  size(A,l); 
per  =  0.95; 
red  =  0.98; 
zstate  =  zstgoal; 
zestim  =  zesgoal; 
zst  =  ones(kappa,l); 
cnt  =  0; 


%  Place  the  discrete  poles  of  the  state  feedback 
% 

while  max(abs(zst))  >  zstgoal  &  cnt<10 
pst  =  (log(zstate)/Ts):... 

(log(per*zstate)Ars  -  log(zstate)/Ts)/(kappa-l):... 
(log(per*zstate)/Ts); 

%pst(10)  =  log(0.1)/Ts 
Frc  =  -place(Ac,Bc,pst); 

Fr  =  -place(A,B,exp(eig(Ac+Bc*Frc).*Ts)); 

zst  =  eig(A+B*Fr); 

if  max(abs(zst))  >  zstgoal 

cnt  =  cnt  +  1 ; 

zstate  =  red  *  zstate; 

else 

cnt=  12; 

end 

end 

pes  =  ones(kappa,l); 
cnt  =  0; 
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% 

%  Place  the  discrete  poles  of  the  estimator 
% 

while  max(abs(pes))  >  zesgoal  &  cnt<10 
pes  =  (log(zestim)Ars):... 

(log(per*zestim)/Ts  -  log(zestim)/Ts)/(kappa-l):... 
(logCper*  zestim)/Ts) ; 

%  pes(lO)  =  log(0.5)/Ts 
Lc  =  -(place(Ac’,Cc’,pes))’; 

L  =  -(place(A,C’,exp(eig(Ac+Lc*Cc).*Ts)))’; 

zes  =  eig(A+L*C); 

if  max(abs(zes))  >  zesgoal 

cnt  =  cnt+  1; 

zestim  =  red  *  zestim; 

else 

cnt  =12; 
end 
end 
% 

%  Place  the  continuous  poles  of  the  state  feedback 
% 

polel  =  -2; 
pole2  =  -3; 

conpoles  =  polel  :(pole2-polel)/(kappa-l):pole2; 

K  =  -place(Ac, Be, conpoles); 
disp(’ Eigen  Values  of  A  +  B*Fr’); 
eig(A+B*Fr) 

disp(’ Eigen  Values  of  A  +  L*C’); 
eig(A+L*C) 

dispC’Eigen  Values  of  Ac  +  Bc*K’); 
eig(Ac+Bc*K) 
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A.6  Optimize.!!! 

function  [v]  =  optimize(cin); 

% 

%  function  [v]  =  optiniize(cin); 

% 

%  This  m-file  will  run  a  constrained  quadratic  optimization  to  find  the 
%  quasi-reference  v(k)  that  will  be  input  to  the  plant. 

%  INPUTS 

% 

%  cin  -  input  vector  that  is  multiplexed  by  simulink  file 
%  cin(l)  =  t,  time 
%  cin(2)  =  Ts,  time  step 
%  cin(3)  =  kappa.  Number  of  states 
%  cin(4)  =  xi,  Number  of  control  inputs 
%  cin(5)  =  eta,  Number  of  outputs 

%  cin(6)  =  kapfl6,  Number  of  subfl6  states  for  F-16  Model 
%  cin(7)  =  p,  State  Prediction  Horizon 
%  cin(8)  =  q.  Control  Horizon 
%  cin(9)  =  r.  Optimization  Horizon 
%  cin(10:xi+9)  =  u,  control  inputs 
%  cin(10+xi:xi+kapfl6+9)  =  xO,  Initial  F-16  Model  states 
% 

%  Written  by  Chris  Shearer,  Aug  97 
% 

t  =  cin(l)  %  Time 

Ts  =  cin(2);  %  Time  Step 

kappa  =  cin(3);  %  Number  of  Estimator  States 

xi  =  cin(4);  %  Number  of  Control  Inputs 

eta  =  cin(5);  %  Number  of  System  Outputs 

kapfl6  =  cin(6);  %  Number  of  F-16  Model  States 

p  =  cin(7);  %  State  Prediction  Horizon 

q  =  cin(8);  %  Control  Horizon 

r  =  cin(9);  %  Optimization  Horizon 

rho  =  max([p  q]); 

% 

%  Pick  apart  the  optimizer  input  from  Simulink 
%  based  upon  the  order  of  the  multiplexier  (Mux) 

% 

xhat  =  cin(10:kappa+9); 

y  =  cin(kappa+10:kappa+eta+9); 

u  =  cin(kappa+eta+10:kappa+eta+xi+9); 

% 
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%  Now  based  upon  the  current  value  of  y  go  get 
%  A,  B,  C,  D  Discretized  State  Space  Model 
%  Fr,  State  Gains 
%  L,  Estimator  Gains 

%  Zr,  Optimized  v  input  gains  (typically  Identity) 

%  T  and  S,  matrices  used  in  qp  probelm 

%  cD  and  cE,  contraint  equations  for  qp  problem 

% 

[A3,C,D,Fr,L,Zr,T,S,cD,cE,l,m,n,ueq]  =  matget(y); 

% 

%  Now  determine  the  current  trajectory  based  upon 
%  the  current  time,  prediction  and  optimization  horizions 
% 

[s,sinf]  =  trajpnt(t,Ts,p,r,0); 

% 

%  This  is  for  addaptive  control 
% 

%  First  I  get  absolute  max  and  minimum  control  inputs 
%  Then  apply  scalling  equations  (ail  =  -0.0600  *  ail/rollrate) 
%  add  a  1.25%  buffer 
%  then  apply  new  max  an  mins 
% 

if  1  ==  2 
bank  =  y(l); 

bankdotp  =  (s(l-i-eta)-s(l))/Ts; 

bankdote  =  (s(l-i-eta)-bank)/T^; 

altdot  =  (s(2+eta)  -  s(2))/Ts; 

altddotp  =  (s(2+2*eta)  -  2*s(2+eta)  +  s(2))/Ts'^2; 

alt  =  y(2); 

altddote  =  (s(2+2*eta)  -  2*s(2+eta)  +  alt)/Ts'^2; 

Umax  =  m(5:8)’  +  ueq; 
if  bankdotp  ==  0  &  altdot  ==  0 
uail  =  (-0.06/2)*bankdote; 
buffa  =  (0.5/1 00)*umax(3); 
elseif  altdot  ==  0 
uail  =  -0.06*bankdotp; 
buffa  =  (0.5/100)*umax(3); 
else 

uail  =  -0.06*bankdotp  -i-  (-0.06/2)*bankdote 

buffa  =  (0.5/100)*umax(3); 

end 
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uailmax  =  uail  +  buffa; 
uailmin  =  uail  >  buffa; 
load  ail 

ailu  =  [ailu  uailmax]; 
aill  =  [aill  uailmin]; 
save  ail  ailu  aill 
m(3)  =  -uailmin; 
m(7)  =  uailmax; 
end 

if  1  ==  2 

uell  =  -altddotp/45; 
uel  =  uell; 

uelmax  =  uel  +  (2/100)*umax(2); 
uelmin  =  uel  -  (2/100)*umax(2); 
m(2)  =  -uelmin; 
m(6)  =  uelmax; 
end 

[vinf]  =  vinfy(A,B,C,Fr,Zr,sinf,r,rho); 

'Rtar  =  [  xhat’  y’  vinf’  s’  ]*T, 

[x, lambda, how]  =  qp(2*S,Tstar,cD,cE*[xhat’  vinf’  y’  u’  1’  m’  n’]’); 

%  Unconstrained 

%[x, lambda, how]  =  qp(2*S,l^tar,0*cD,cE*[xhat’  vinf’  y’  u’  1’  m’  n’]’ 
v  =  x(l:xi); 
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APPENDIX  B  -  Operating  Envelope  and  F-16  Planform  and 

Longitudinal  Views 

B.l  F-16  Planform  and  Longitudinal  Views 
Removed  for  Distribution  Purposes 
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B.2  F-16  Operating  Envelope/Tiirn  Performance  -  Sea  Level 


Removed  for  Distribution  Purposes 
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APPENDIX  C  -  LTI  Models 


C.1  HARV 


A 


B 


C 
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and  D  matrices  that 
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27 


were  used  in  the  simulations  were  scaled  using  the  following 


relationship 

Asci  =  T-^AT 
Bsci  =  T-^AS 
Csci  =  R-^AT 
Dsci  =  R~^AS 
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C.2  F-16 


Listed  below  are  the  LTI  state  space  matrices  used  in  the  linear  simulations 


-2.53e  -  02 

3.94e  +  01 

-3.22e  +  01 

2.41e  -  01 

4.28e  -  01 
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APPENDIX  D  -  Extra  Simulation  Results 


D.l  Combined  35%  C.G.  Location,  Linear  vs  Nonlinear  Comparison 


Figure  94.  Linear  and  Nonlinear  Comparison  Figure  95.  Linear  and  Nonlinear  Comparison 
-  Longitudinal  Control  Inputs  -  Lateral  Control  Inputs 


C.G.  =  0.35  C.G.  =0.35 


Figure  96.  Linear  and  Nonlinear  Comparison 
-  Bank  Angle 


Figure  97.  Linear  and  Nonlinear  Comparison 
-  Side  Slip  Angle 
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Velocity  (fVsec) 


Figure  98.  Linear  and  Nonlinear  Comparison  -  Pitch  Rate 


Figure  99.  Linear  and  Nonlinear  Comparison  Figure  100.  Linear  and  Nonlinear  Compari- 
-  \felocity  son  -  Altitude 


Figure  101.  Case  3  -  Bank  Angle  Output 
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Figure  103.  Case  3  -  Longitudinal  Controls 


Longitudinal 


Figure  128.  Case  19  -  Longitudinal  Control  Figure  129.  Case  19  -  Lateral  Control  Inputs 
Inputs 


D.9  Case  21 


Figure  130.  Case  21  -  Altitude  Output 


Figure  131.  Case  21  -  \felocity  and  Bank  An¬ 
gle  Outputs 
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Figure  132.  Case  21  -  Beta  and  Theta  Outputs 
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Figure  133.  Case  21  -  Longitudinal  Inputs 


Figure  134.  Case  21  -  Lateral  Inputs 


D.IO  Case  22 


Figure  135.  Case  22  -  Altit 


Figure  13 


Outputs  Longitudinal  Control  Inputs 


Tims  (sees)  Time  (secs) 

Figure  138.  Case  22  -  Longitudinal  Inputs  Figure  139.  Case  22  -  Lateral  Inputs 

Case  24 


Time  (secs)  Time  (secs) 


Figure  140.  Case  24  -  Altitude  Output 


Figure  141.  Case  24  -  Bank  Angle  Output 


Longitudinal  Control  Inputs 


Time  (secs) 


Figure  142.  Case  24  -  \felocity,  Beta  and  Theta  Outputs 


Figure  143.  Case  24  -  Longitudinal  Controls 


Figure  144.  Case  24  -  Lateral  Controls 


Outputs 


D.12  Case  26 


Figure  145.  Case  26  -  Altitude  Output  Figure  146.  Case  26  -  Bank  Angle  Output 


Figure  147.  Case  26  -  \felocity.  Beta  and  Theta  Outputs 


Longitucfinal  Contrtrf  Inputs 


Figure  148.  Case  26  -  Longitudinal  Controls 


Figure  149.  Case  26  -  Lateral  Controls 
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